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ABSTRACT 


As part of the SEASAT Satellite program of NASA, a 
survey of representative primitive equation atmospheric 
prediction models that exist m the world today was written 
for the Jet Propulsion Laboratory under Contract No. 954547 
(Subcontract under NASA Contract No. NAS7-100, Task Order 
No. RD-151) . Seventeen models developed by eleven different 
operational and research centers throughout the world are 
included m the survey. The surveys are intended to be 
tutorial m nature describing the features of the various 
models in a systematic manner. 
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I. 


Int roduct ion 


At the request of the Jet Propulsion Laboratory and m 
support of the SEASAT satellite program of NASA, Ocean Data 
Systems, Inc. has conducted a survey of representative at- 
mospheric prediction models m use at operational and research 
centers throughout the world. These surveys are intended to 
be tutorial m nature describing the features of the various 
models m a systematic manner. 

When conducting a survey on some subject, one is 
immediately faced with the problem of choosing the items to 
be surveyed. If the survey is to be exhaustive, then the 
only problem is making sure that it is truly exhaustive. 

Since this atmospheric prediction model survey was not in- 
tended to be exhaustive, but to be representative of the 
types of primitive equation atmospheric prediction models 
that exist m the world today, criteria for model selection 
was required. 

This survey has attempted to include models from many 
different organizations, models with a wide range of appli- 
cations - research to operational, models of varying- types 
- gridpomt and spectral, and models of varying areal co- 
verage - global to limited area. In addition, the availa- 
bility of literature describing the models, a desire to avoid 



repetition that would be caused by surveying quite similar 
models, and an attempt for the survey to be representative 
of atmospheric prediction models in the world today, in- 
fluenced the model selections . 

The accompanying chart represents a crude attempt to 
describe the "genealogy" of the various models surveyed. 

The models are denoted by the organizations that developed 
them and have been divided into two basic categories - 
gndpomt and spectral. The solid lines represent direct 
descendents while the dotted lines represent interchange of 
procedures and ideas. While there are not very many direct 
descendents, there is a tremendous amount of interchange of 
ideas within the modelling community. 

Seventeen different models from the eleven different 
organizations shown on the chart were included m the survey. 
Where multiple models from a single organization were included, 
they are ususally quite similar and differ only m numerical 
techniques or m areal extent (e.g., a limited area model). 

A standard format was devised for the model surveys in 
order that features of the various models can more easily 
be compared. That is, the description of the variables, 
areal coverage, physical and dynamical processes, numerical 
methods, and so forth, are presented m the same order, topic 
by topic for each prediction model. 
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Finally, it should be pointed out that some of the 
model surveys are more detailed than others. This, of 
course, reflects the amount of information that was availa- 
ble concerning the model. Also, all of the model descrip- 
tions represent a snapshot of the various prediction models 
at some point m time defined by a combination of the source 
material available and the knowledge of the author concer- 
ning the particular model. While the availability and 
accuracy of the source material heavily impacts the accu- 
racy of a given model description, the ultimate responsi- 
bility for the accuracy of the model description rests with 
the author of the surveys . 

The following Section, Section II, contains the surveys 
of the models included m this atmospheric prediction model 
survey. Section III discusses the inclusion of SEASAT Satel- 
lite data m the forecast models m relation to the analysis 
problem. Mention is also made of the boundary layer para- 
meterization problem as related to the use of SEASAT Satel- 
lite data. 
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II. Atmospheric Prediction Model Surveys 


This section contains the atmospheric prediction model 
surveys. Each model is presented m a standard format as 
a subsection m order to more easily compare the features of 
the various models. 
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III. Summary 


The previous section contains surveys of seventeen 
different atmospheric prediction models developed by eleven 
different organizations. These models are representative 
of the types of primitive equation atmospheric prediction 
models that exist m the world today. Included are models 
from many different organizations, models with a wide range 
of applications - research to operational, models of varying 
types - gridpomt and spectral, and models of varying areal 
coverage - global to limited area. 

In the descriptions of the models, the initial data is 
either treated as input to the program m the case of the 
forecast models, or as being self-generated m the case of 
the general circulation (climate simulation) models. The 
descriptions of the inputs to the forecast models usually 
make reference to an analysis procedure which produces data 
for use by the forecast model. Thus, the model descriptions 
only very indirectly make reference to the analysis problem. 
This is, of course, reasonable since discussion of analysis 
procedures is beyond the scope of this survey. However, m 
order to make a forecast an initial state is required and 
some mention should be made of the analysis problem including 
data sources. 
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Historically, the analysis data cycle has been geared 
to those times at which upper-air observations (rawinsondes) 
are made (0000 and 1200 GMT, mainly) . Analyses of perti- 
nent atmospheric parameters are made at these times for the 
entire three-dimensional domain of interest. However, with 
the availability m recent years of satellite temperature/ 
moisture profiles (SIRS, VTPR) derived from radiance mea- 
surements, considerable attention has been given to methods 
for incorporating these "off-time" data into atmospheric 
analysis structures. 

The satellite data available today is spread through 
both space and time. Satellites m polar orbits produce 
data m strips corresponding to orbits around the earth and 
the width of the field of view of the sensors. Satellites 

i 

m other types of orbits also produce data that is distri- 
buted m space and time m a manner dependent upon the orbit 
and the sensor characteristics; but, not as the classical 
analysis cycle would prefer it - global coverage at the 
canonical times of 0000 and 1200 GMT. 

SEASAT will add to the amount of data available for 
analysis by producing large amounts of data over the oceans 
including wave heights, surface winds and ocean temperatures. 
This data will, of course, be "off- time" data as is other 
satellite data and methods will have to be found to incor- 
porate the data into the analysis structure. 
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In addition , SEASAT data will all be at the ocean sur- 
face which will introduce a different class of problem m 
the utilization of the data. With the exception of variables 
such as surface pressure, forecast models carry data at 
levels above the surface. In fact, m most models, surface 
values of many variables are obtained through diagnostic 
boundary layer calculations. SEASAT data will present the 
inverse of this problem - many surface values will be avail- 
able from measurements and boundary layer parameters zations 
must be devised to couple this information into the levels 
above the surface. 

The previous paragraphs have given a brief discussion 
of the analysis problem since the forecast problem cannot be 
addressed without at least mentioning the related analysis 
problem. In particular, the problem of incorporating "off- 
time" , strip data and data at "inconvenient" levels into the 
analysis structure has been pointed out. These problems are 
all being actively investigated and the solutions or, at least, 
partial solutions will certainly be included as part of future 
atmospheric analysis/prediction systems. 
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A. AFGWC 6- Layer Hemispheric Operational Model 

The Air Weather Service atmospheric prediction model 
(AWSPE) is currently executed twice daily to produce three 
day forecasts for the Northern Hemisphere. The model is 
an adaptation of the NMC model without moisture (Flattery, 
1975) . Since the forecast model is continually being 
modified, the description that follows is, thus, a snapshot 
of the state of the model as of mid-1976. 

MODEL DESCRIPTION 

1.0 Variables 

1 . 1 Independent 

The independent variables are x and y derived from a 
polar stereographic map projection in the horizontal, sigma 
m the vertical defined as : 



where P is the pressure at a given quasi-horizontal surface 
above some point and P_ is the pressure at a surface below, 

J-i 

and time, t. 
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1. 2 Dependent 


1.2.1 Prognostic 

The prognostic variables are the x wind component, u; ' 
the y wind component, v; the potential temperature, 8; and 
the vertical pressure gradient, — , at the three interfaces 
(between the planetary boundary layer and the troposphere, 
between the troposphere and the stratosphere, and between 
the stratosphere and the constant 0 layer) as shown m 
Figure A-2. 

1.2.2 Diagnostic 

The diagnostic variables are the geopotential , $ , 
and the vertical velocity, w = -a. 

2.0 Domain 

2.1 Horizontal 


The horizontal domain is a polar stereographic projec- 
tion (true at 60 °N) of the Northern Hemisphere using x-y 
coordinates with map factor 


1 t sin 60° 
1 + sin. (|> 


[A. 2] 


The grid resolution is '381 km at 60 °N and the grid 
dimensions are 53 x 57 which covers virtually the entire 



Northern Hemisphere as shown m Figure A-l. (The standard 
octagon is also shown.) There is no staggering of variables 
m the horizontal. 

2.2. Vertical 

The vertical coordinate is sigma as defined by Equation 
[A.l]. Figure A-2 shows the vertical structure of the 
model with the four distinct sigma domains. The variables 
u, v and 0 are carried at the midpoints of the various 
sigma layers and the variables $ , 77 and p are carried at 
the layer interfaces. The vertical velocity, a, is defined 
at layer interfaces, but at the center of a horizontal 
gridpoint square.. 
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AWS Northern Hemisphere Polar Stereographic Grid 
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Figure A- 2 

Vertical Structure of the AWS Atmospheric Prediction Model 
2 . 3 Temporal 

■» 

Two time levels of data are required (except during 
the startup step) since a centered time differencing scheme 
is used with a twenty minute time step. There is no 
staggering of variables m time. 
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3 . 0 Governing Equations 


The governing partial differential equations, written 
for an x-y map projection with map factor m, are listed 
below. 

Wind component m the x direction: 
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Wind component m the y direction: 
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Thermodynamic energy equation: 


30_ 

at 


= -m [ u 


36 

ax 


+ v 


30 . i 
3y J 


aa. + 

aa 


H 


Pressure tendency equation: 
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Hydrostatic equation: 


= _c 0 3-B- 
3a P 3a 


Exner function: 


ir 


( 


1000 


) 


K 


In the above equations , 


u = x direction wind component 

v = y direction wind component 

0 = potential temperature 

$ = gz = geopotential 

p = pressure 

f = Coriolis parameter 

a = = vertical velocity 

Fx,Fy = stress components 


H 

C 

P 

R 

K 


= diabatic heatmg/coolmg term 

= specific heat of dry air at constant pressure 

= gas constant for dry air 
R 


[A. 7] 


[A. 8] 
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4.0 Numerics 


Finite difference methods are used to approximate 
the solution of the forecast model equations . 

4.1 Spacial Differencing 

Centered, second-order differencing using Shuman's 
(1968) horizontal space averaging technique is used to 
approximate the x and y derivatives m the various sigma 
domains. The vertical derivatives are also approximated 
by second-order differences. 

4.2 Temporal Differencing 

Centered time differencing with Shuman's pressure 
gradient force time averaging is used to obtain a time 
step of twenty minutes. The computation is started with 
a forward time step. 

4 . 3 Boundary Conditions 

The horizontal boundary conditions consist of insula- 
ted, slippery walls placed on the grid as shown m 
Figure A-l. That is, no heat or momentum flow is 
permitted through the wall^ but a parallel flow rs per- 
mitted. along the wall.. 
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The vertical boundary condition that a = o 
for a material surface is specified: at the top of the 

model; the boundary surface between the topmost compu- 
tational layer and the stratosphere; the tropopause; 
and the earth ' s surface . As a consequence of the above 
approach, no mass exchange takes place among the isen- 
tropic, stratospheric and tropospheric domains. However, 
mass exchange can take place between the boundary layer 
and the troposphere. 

5.0 Processes and Effects 

5.1 Dynamics 

Terms representing the horizontal and vertical 
advection and convergence/divergence of heat and momentum 
are included as are terms representing the pressure 
gradient and Coriolis forces m the momentum equations . 

Orography that is rough but consistent with the 
grid resolution is included to add realism to the flows. 

Horizontal stress according to Shuman (1968) is 
applied at the lowest level to simulate vertical momentum 
transport at the surface. 
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5.2 Heat Sources/Sinks 


Sensible heating m the planetary boundary layer 
over the oceans is included. The ocean temperature is 
spacially varying, but held constant during the forecast. 

A dry convective adjustment procedure is included 
to preclude hydrostatic instability of the model 
atmosphere . 

5.3 Moisture Sources/Sinks 

There are none since a dry, adiabatic fluid is 
considered. 

6 . 0 Inputs 

The forecast model initial data is obtained from 
hemispheric analyses of the mass structure (temperature 
and geopotential) between 1000 mb and 100 mb. The 
analyses are performed by a program called MULTAN which 
is the AFGWC objective analysis program and is based on 
the work of NMC. 

The initial wind fields are derived from the linear 
balance equation solved on the mandatory pressure surfaces. 
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The analyzed and derived fields on pressure 
surfaces are then interpolated (linearly m lnp) to 
obtain initial data on the model sigma surfaces. 

The sea surface temperature is obtained from Fleet 
Numerical Weather Central analyses. 

7 . 0 Outputs 

Fields on the standard pressure surfaces are produced 
by interpolation from the model sigma surfaces to the 
pressure surfaces during an output cycle. Various filters 
are used to render the output fields more visually 
pleasing. 

8.0 Computer Information 

The forecast model is run on a Umvac 1110 Computer 
System and is normally executed twice daily. The model 
is largely programmed m Fortran V (Umvac dialect of 
FTN IV) with some routines in assembly language. 

9 . 0 References 

Flattery, T.W.,"The Air Weather Service Primitive Equation 
Model", AFGWC User Documentation, May 1975. 

Shuman, F.G. and J.B. Hovermale, "An Operational Six-Layer 
Primitive Equation Model" , Journal of Applied Meteor- 
ology, Vol. 7, No. 4, August 19&8, pp. 525-547. 
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B. FNWC 5-Layer Hemispheric Operational Model 

The FNWC operational atmospheric prediction model is 
executed twice daily to produce three day forecasts for 
the Northern Hemisphere. The model was designed and 
developed for operational use by Kesel and Winnmghoff 
(1972) during the period 1969-72. Since that time many 
changes have been made to the model to improve the quality 
of the forecasts produced and to better utilize the 
available computer resources. The description that follows 
is, thus, a snapshot of the state of the model as of mid- 
1976. 

MODEL DESCRIPTION 

1.0 Variables 

1 . 1 Independent 

The independent variables are x and y derived from a 
polar stereographic map pro 3 ection in the horizontal, sigma 
m the vertical defined as 

a = [B.l] 

IT 

where it is the local terrain pressure, and time, t. 
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1 . 2 Dependent 


1.2.1 Prognostic 

The prognostic variables are the terrain pressure, ir , 
the temperature, T, the water vapor pressure, q, the x 
wind component, u, and the y wind component, v. A secondary 
prognostic variable is the precipitation, Pr, obtained 
from both the large scale precipitation algorithm and the 
cumulus parameterization. 

1.2.2 Diagnostic 

The diagnostic variables are the geopotential r $ , 
and the vertical velocity, w = -o . 

2 . 0 Domain 


2.1 Horizontal 


The horizontal domain is a polar stereographic projection 
(true at 60 °N) of the Northern Hemisphere using x-y coordi- 
nates with map factor: 


l+sm6Q ° 
l+sin<J> 


[B . 2 ] 


The horizontal grid resolution is 381 km at 60°N resulting 
m a 63 x 63 grid as shown m Figure B-l. There is no 
staggering of variables m the horizontal. 
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Figure B-I: Northern Hemisphere Polar Stereographic Grid. 

The pole to equator distance is 31.205 grid lengths. 
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2.2 Vertical 


The vertical coordinate is sigma as defined by Equation 
[B.l], Figure B-2 shows the model sigma surfaces and the 
placement of the variables in the vertical. 
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Figure B-2: Vertical Structure of the FNWC 

5-layer Operational Forecast Model . 


2 . 3 Temporal 

Two time levels of data are required (except during 
the startup step) since a centered time differencing scheme 
is used with a ten minute time step. There is no staggering 
of variables m time. 
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3.0 Governing Equations 

The governing partial differential equations, written in 
flux form for an x-y map projection with map factor m, are 
listed below. 

Momentum equation m the x direction: 


0 (iru) 2 r 3 ,uutk , 3 ,uvir. . 

~Tt = ~ m <— > + w ( -ir' ] 

+ uvf - m (ttJ^ + RT + F 
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Momentum equation m the y direction: 
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Thermodynamic energy equation: 
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Moisture conservation equation: 
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Pressure tendency equation: 
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Hydrostatic equation: 

9$ RT 

3o o 


[B . 8] 


In the above equations , 

u = x direction wind component 

v = y direction wind component 

it = terrain pressure 

T = temperature 

q = vapor pressure 

w = _ vertical velocity 

f = Coriolis parameter 

F ,F = stress components 

x y 

D u ,D v ,D T ,D^ = diffusion terms 
H = diabatic heating/cooling term 

Q = moisture source/sink term 

# = gz = geopotential 

R = gas constant for dry air 

Cp = specific heat of dry air at constant pressure 
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4.0 Numerics 


Finite difference methods are used to approximate the 
solution of the forecast model equations. 

4.1 Spacial Differencing 

Centered, fourth-order differencing of the horizontal 
advection and pressure gradient force terms is used to 
reduce the phase propagation error. All other derivatives 
are approximated by centered, second-order differences. 

The lateral diffusion terms shown m the momentum and 
thermodynamic energy equations take the form of diffusion of 
the departures from the initial state, while the lateral 
diffusion term m the moisture equation is based on the 
diffusion of the moisture only. Thus, the term D T is 
defined as : 

D = ^ 2 A 2 (tt°T° - tt t t x ) [B.9] 

i a 

2 

where K is the diffusion coefficient, A represents the Lap- 
lacian operator, the superscript zero denotes the initial 
value and the superscript x denotes the current value. 

4.2 Temporal Differencing 

Centered time differencing with Shuman's pressure- 
gradient force time averaging is used to offset the more 
stringent linear stability requirements of fourth- order 
differencing resulting in a time step of ten minutes. 
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Robert time filtering of the solutions is also done to 
control computational noise. 

The computation is started with an Euler-backward (or 
Matsuno) step. This procedure is repeated every six hours 
after each output cycle. 

4.3 Boundary Conditions 

The horizontal boundary conditions are of the restora- 
tion type; that is, persistence below 7.5°N, a blend of 
persistence and prediction between 7.5° and 15 °N, and total 
prediction above 15 °N. 

The vertical boundary conditions are w = 0 at a = 0 
and 1 . 

5.0 Processes and Effects 

5 . 1 Dynamics 

Terms representing the three-dimensional flux 
divergence of heat, moisture and momentum are included 
as are terms representing the pressure gradient and 
Coriolis forces m the momentum equations. Adiabatic 
temperature changes are also modeled. 

The Kunhara local pressure surface technique is used 
to calculate the pressure gradient force terms m order to 
reduce computational noise over high, irregular terrain. 

Orography that is rough but consistent with the grid 
resolution is included to add realism to the flows. The 
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terrain gradients are limited to 2000 meters per grid 
interval and the terrain is smoothed. 

Friction according to Shuman is included in the lowest 
layer (o = 1.0 to 0.8). 

5.2 Heat Sources/Sinks 

Both solar and terrestrial radiation are included m 
a manner similar to that of the Mmtz and Arakawa model of 
the late sixties. Moisture is included only m the lower 
three sigma layers, and heating rates are evaluated for 
two gross layers only, from o = 1.0 to 0.6 and from a = 0.6 
to 0.2. Clouds are included at the a - 0.7 level using a 
parameterization due to Smagonnsky. 

Sensible heating m the Arakawa type constant flux, 
gross boundary layer is evaluated for the lowest sigma 
layer and a land-sea-ice discriminator is included. The 
ocean temperature is spacially varying, but held constant 
during the forecast. 

Release of latent heat due to large scale precipita- 
tion is included m the moisture carrying layers. 

A parameterization of cumulus ensembles (moist convec- 
tive adjustment) based on the 1968 version by Arakawa et alL 
redistributes heat and moisture m the vertical when the 
appropriate conditions prevail (conditional instability 
of the layers) . A dry convective adjustment procedure 
is included to preclude hydrostatic instability. 
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5.3 Moisture Sources/Smks 


The large scale precipitation process removes super- 
saturated states from the model atmosphere and may result 
m evaporation m layers below the supersaturated layer. 

The moist convective adjustment process can result m 
the evaporation/condensation of moisture m the lower layers. 

Evaporation from the surface in the constant flux, 
gross boundary layer is also included m the lowest sigma 
layer. 

6.0 Inputs 

The forecast model initial data is obtained from hemis- 
pheric analyses of the mass structure between 1000 mb and 
100 mb. The analyses are performed using a technique called 
Fields by Information Blending (FIB) devised by Holl (1972) . 
The 50 mb level heights and temperatures are obtained through 
linear regression formulae using the 100 mb heights and tem- 
peratures. Dewpoint depression, sea level pressure and sea 
surface temperature analyses are also made using the FIB 
technique. 

The rotational component of the winds is derived from 
the solution of the balance equation on pressure surfaces m 
which the nonlinear terms are approximated geostrophically . 
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The analyzed fields on pressure surfaces are inter- 
polated (linearly in £n(p)) to obtain initial data on the 
model sigma surfaces. 

7 . 0 Outputs 

A normal output cycle is made every six forecast hours. 
The sea level pressure, precipitation, and heights and 
temperatures at standard pressure levels from 1000 mb to 
100 mb are output. The fields on pressure surfaces are 
produced by interpolation (or extrapolation algorithms 
below the terrain) from the sigma surface data. Various 
filters are used to render the output fields more visually 
pleasing . 

Numerous other fields (winds, heating rates, etc) are 
available during the output cycle. 

8 . 0 Computer Information 

The forecast model is run on two Control Data Corporation 
6500 computer systems (four processors) coupled through 
extended core memory. (See Morenoff, et al. 1972.) The 
model is normally executed twice per day and requires 
about ninety minutes to produce a seventy-two hour forecast. 

The software system is slightly modified CDC and the 
model is written mainly m Fortran with some routines 
m assembly language. Various special purpose routines 
peculiar to FNWC are also used by the model. 
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C. FNWC 5-Layer Global Model 

The FNWC global primitive equation atmospheric 

prediction model has been proposed as a successor to the 

* 

current hemispheric operational model. (See Section B.) 
Development of this model was begun by Winninghoff (1971) 
and has been continued by Elias (1973) , Mihok (1974) and 
Mihok and Kaitala (1976) . Since this global model is an 
outgrowth of the operational hemispheric model, there 
are many similarities between the two models. The grid 
structure and finite difference techniques used are quite 
different, while most of the physical parameterizations 
other than the treatment of the boundary layer are 
identical . 

MODEL DESCRIPTION 

1 . 0 Variables 

1.1 Independent 

\ 

The independent variables are latitude, 6, and longitude, 
A , of a spherical polar coordinate system in the horizontal, 
sigma m the vertical defined as 

cr = 2. [c.l] 

TT 

where ir is the local terrain pressure, and time, t. 
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1.2 Dependent 


1.2.1 Prognostic 

The prognostic variables are the terrain pressure, it, the 
temperature, T, the water vapor pressure, q, the latitudinal 
wind component, u, and the longitudinal wind component, v. 

A secondary prognostic variable is the precipitation, Pr , 
obtained from both the large scale precipitation algorithm 
and the cumulus parameterization. 

1.2.2 Diagnostic 

The diagnostic variables are the geopotential , 
and the vertical velocity, w = -cr. 

2 . 0 Domain 

2.1 Horizontal 

The horizontal domain is a global latitude, longitude 
grid with 0 and A as the coordinates. The resolution is 
5° between values of like variables giving a 71 x 72 grid 
with a grid distance of 556 km at the equator. The poles 
are considered as mass points and are treated separately-. 

The variables are staggered in the horizontal as shown m 
Figure C-l with velocity components of the wind carried at 
one set of points- and mass variables at the other. 
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Figure C-l 


Horizontal Grid Structure of the FNWC 5-Layer Global Model 


e-3 






2 . 2 Vertical 


The vertical coordinate is sigma as defined by Equation 
[C.13. Figure C-2 shows the model sigma surfaces and the 
placement of the variables m the vertical. 
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Vertical Structure of the FNWC 5-Layer Global Model 
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2 . 3 Temporal 

Two time levels of data are required (except during 
the first time step) since a centered time differencing 
scheme is used with Shuman's pressure gradient force aver- 
aging technique and Robert time filtering. A nine minute 
time step is used. There is no staggering of variables in 
time. 
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3 . 0 Equations 


The governing partial differential equations written 
m flux form for a spherical polar coordinate system, 0 
arid A, m the horizontal and a sigma coordinate m the ver- 
tical are listed below. 

Momentum equation m the latitudinal direction: 


3 (ttu) 
9 1 


acos0 


[ 


3 (uutt) . 3Juvrrcos0) 


3A 


30 
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3 (wu) 
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+ (f + 


utan0 


) ttv - 
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tt3$ 
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+ F + D 
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Momentum equation in the longitudinal direction: 


3 (irv) 
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Thermodynamic energy equation: 


9 (ttT) 
3t 
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Moisture conservation equation: 

3 ( TTq) _ _ 1 r 3 ( 7Tuq) , 3 (7rvqcos6) ■■ 3 (wq) 

3t acos0 1 3A 30 3a 

+ Q-ft + D q 

Pressure tendency equation : 

3 it f 1 1 r 3(uir) r 3(vTrcos0)., 

3t “ J 0 acos0 L 3A 30 J ao 

Hydrostatic equation: 


3$ 

RT 

3a 

a 

the above 

equations : 

u = 

zonal wind component 

v = 

meridional wind component 

7T = 

terrain pressure 
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[C. 4] 

EC. 5] 
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R = gas constant for dry air 

C = specific heat of dry air at constant pressure 

hr 

q = vapor pressure 

w = - ^ = vertical velocity 

a = radius of the earth 

f = Coriolis parameter 

F x ,Fy = skin friction 

D^ = lateral diffusion of momentum 
m 

D t = lateral diffusion of heat 
Dg = lateral diffusion of moisture 

H = diabatic heating/cooling term 

Q = moisture source/sink term 

$ = gz = geopotential 

4 . 0 Numerics 

Finite difference methods are used to approximate the 
solution of the forecast equations. 

4 . 1 Spacial Differencing 

Centered, fourth-order differencing in the horizontal 
is used to more accurately compute the derivatives and 
reduce the phase error. Vertical derivatives are 
approximated by centered, second-order differences. Due 
to the staggering o'f variables m the horizontal, space 
averaging of variables is sometimes required in the 
finite difference approximations. 'This renders the difference 
scheme slightly dissipative and, hence, somewhat more stable 
than a- non- space averaged scheme. . 
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The polar calculations follow Arakawa m that the change 
of the variables defined at the poles (ir r T and q) is a re- 
sult of the meridional mass flux at all points on the latitude 
circle around the particular pole. 

The lateral diffusion terms shown in the momentum, thermo- 
dynamic energy and moisture equations take the form of dif- 
fusion of the departures from the initial state. Thus, the 
term D T is defined as: 

D = - j A 2 (tt°T° - tt X T t ) [C . 8] 

acosSAx^ 

2 

where K is the diffusion coefficient, A represents the 
Laplacian operator, the superscript zero denotes the initial 
value and the superscript t denotes the current time. 

4.2 Temporal Differencing 

Centered time differencing with Shuman's pressure- 
gradient force time averaging is used to allow a longer time 
step of nine minutes to be used. This procedure is coupled 
with Robert time filtering of the solutions to help control 
computational noise. 

The computation is started with an Euler-backward or 
Matsuno step. 

The removal of the overspecification of the tendencies, 
caused by the decreasing grid size as one approaches the 
poles, is accomplished through the use of a modified Arakawa 
three point averaging scheme. 
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4 . 3 Boundary Conditions 


The horizontal boundary conditions are periodicity m 
longitude and the poles are treated as mass points as ex- 
plained m Section 4.1. 

The vertical boundary conditions are w = 0 at o=Q and 1. 

5.0 Processes and Effects 

5 . 1 Dynamics 

Terms representing the horizontal and vertical advection 
and convergence/divergence of heat, moisture and momentum are 
included as are terms representing the pressure gradient and 
Coriolis forces m the momentum equations. Adiabatic tem- 
perature changes are also modeled. 

Orography that is consist ent with reality and the grid 
resolution is included to add realism to the flows . The 
terrain gradients are limited and the terrain is smoothed. 

The skin friction, F and F , is included m the lowest 

x y 

sigma layer and is calculated from the planetary boundary 
layer (PBL) parameterization devised by Kaitala. This PBL 
parameterization is based on the premise that the steady 
state, barotropic PBL is completely determined by a pair of 
external parameters: the surface Rossby number, and the 

stability parameter. 
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5 . 2 Heat Sources/Sinks 


Both solar and terrestrial radiation are included in a 
manner similar to that of the Mmtz and Arakawa model of the 
late sixties. Moisture is included only m the lower three 
sigma layers, and heating rates are evaluated in two gross 
layers only, from cr = l.o to 0.6 and from a = 0.6 to 0.2. 
Clouds are included at the a = 0.7 level using a parame- 
terization due to Smagonnsky. 

Sensible heating m the PBL is evaluated for the lowest 
sigma layer, and a land-sea-ice discriminator is included. 

The ocean temperature is spacially varying, but held constant 
during the forecast. 

Release of latent heat due to large scale precipitation 
is included m the moisture carrying layers. 

A parameterization of cumulus ensembles (moist convec- 
tive adjustment) based on the Arakawa method and modified 
somewhat by Kaitala redistributes heat and moisture m the 
vertical when the appropriate conditions prevail. A dry 
convective adjustment procedure is included to preclude 
hydrostatic instability. 

5.3 Moisture Sources/Sinks 

The large scale precipitation process removes super- 
saturated states from the model atmosphere and may result 
in evaporation m layers below the supersaturated layer. 
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The moist convective adjustment process can result m 
the evaporation/condensation of moisture m the lower layers. 

Evaporation from the surface is evaluated as part of the 
PBL computations and included in the lowest sigma layer. 

6 . 0 Inputs 

The forecast model initial data is obtained from a global 
2 1/2° objective analysis of the mass structure between 1000 mb 
and 100 mb. The analyses are performed using a technique 
called Fields by Information Blending (FIB) devised by Holl 
(1972). The 50 mb heights and temperatures are obtained 
through linear regression formulae from the 100 mb heights 
and temperatures. Dewpoint depression, sea level pressure 
and sea surface temperature analyses are also made using the 
FIB technique. 

Winds (rotational component) are calculated on pressure 
surfaces from the stream function (i^) obtained from the solu- 
ton of the full balance equation m which the nonlinear terms 
are approximated geostrophically . The balance equation is 
solved first m the southern hemisphere while holding the 
pole and the equator plus one latitude circle to the north 
constant. This extra latitude circle has the same sign and 
magnitude of the Coriolis force as that of the latitude circle 
to the south. The process is then reversed for the northern 
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hemisphere. The winds at the equator are averages of the 
components at the latitude circles adjacent to the equator. 

The fields on pressure surfaces are then interpolated 
(linearly m £n(p)) to the model sigma surfaces. 

7 . 0 Outputs 

A normal output cycle is made every six forecast hours. 
The sea level pressure, precipitation and heights and tem- 
peratures at standard pressure levels from 1000 mb to 100 mb 
are output. The fields on pressure surfaces are produced by 
interpolation (or extrapolation algorithms below the terrain) 
from the sigma surface data. Various filters are used to 
render the output fields more visually pleasing. 

8 . 0 Computer Information' 

The forecast model is programmed mostly m CDC Fortran 
Extended (a dialect of FTN IV) with a few routines m assem- 
bly language. The current computer systems used are a 
CDC 6500 and a CDC Cyber 175. 
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D. GFDL 9 -Layer Global Model 


This Geophysical Fluid Dynamics Laboratory (GFDL) 
atmospheric general circulation model was developed for 
use m long range prediction and general circulation 
simulation experiments. The model is based on earlier 
versions developed at GFDL, but uses a different finite 
difference grid and approximations and contains improved 
par ameter i z ations of the subgrid scale processes. What 
follows is a description of the model as of 1976 as 
described m GARP publication No. 14. 


MODEL DESCRIPTION 

1.0 Variables 

1 . 1 Independent 

The independent variables are the longitude, X, and the 
latitude, <i>; sigma defined as 

a = E [D . 1] 

3 ?*. 

where p* is the surface pressure; and time, t. 
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1 . 2 Dependent 


1.2.1 Prognostic 

The prognostic variables are the longitudinal 
component, u, and the meridional component, v, of the 
wind; the temperature, T; the surface pressure, p*; and 
the moisture mixing ratio, r. Secondary prognostic 
variables are the surface temperature over land, T g/ 
and the soil moisture, W. 

1.2.2 Diagnostic 

The diagnostic variables are 

* 

the vertical cr velocity, a = — , 
velocity, a). 

2 . 0 Domain 

2.1 Horizontal 

The horizontal domain is a global latitude-longitude 
grid with A and as coordinates. The variables are 
specified at the center of the latitude- longitude boxes 
and, thus, the poles are not grid points. There is no 
staggering of variables m the horizontal. Different 
versions of the model exist with varying horizontal 
resolution ranging from a low resolution version with 


the geopotential , 
and the vertical p 
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nineteen points between the equator and the poles (N19) 
and 64 points around a latitude circle to much finer 
resolution versions. 


2 . 2 Vertical 

The vertical coordinate is sigma as defined by 
Equation [D.l]. Figure D-l shows the vertical structure 
of the model and the placement of variables m the vertical . 
Table D-l gives the distribution of the sigma levels and 
their approximate height above sea level. 
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Figure D-l 

Vertical Structure of the GFDL 9-Layer Global Model 
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Table D-l 


Location of the a Levels m the GFDL 9-Layer Model 
and Their Approximate Elevations Above the Surface 
Points at Sea Level, (p* is assumed to be 1000 mb.) 


Level 

<^P/P* 

Height (m) 

0.5 

0.00000000 

CO 

1.0 

.01594441 

27,900 

1.5 . 

.04334139 

21,370 

2.0 

.07000000 

18,330 

2 5 

.11305591 

15,290 

3.0 

.16500000 

12,890 

3.5 

.24081005 

10,500 

4.0 

.31500000 

8,680 

4.5 

.41204675 

6,860 

5.0 

.50000000 

5,430 

5.5 

,60672726 

4,010 

6.0 

.68500000 

3,060 

6.5 

.77337055 

2,110 

/.0 

.83500000 

1,490 

7.5 

.90154066 

860 

8.0 

.94000000 

520 

8.5 

.98010000 

170 

9.0 

.99000000 

80 

9.5 

1.00000000 

0 


2 . 3 Temporal 


Two time levels of data are required (except during the 
startup step) since a centered time differencing scheme is 
used with a ten minute time step for the N19 version of 
the model. There is no staggering of variables m time. 
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3 . 0 Governing Equations 


The governing partial differential equations written 
m flux form and spherical polar coordinates with a 
vertical coordinate sigma are listed below: 

Wind component m the latitudinal direction: 


9t 


9 (p*uu) 
a cos 4>3A 


3 (p*cos<J>uv) 3{crp*u) 


a cos <pd(f> 


3 a 


a + u tan h 
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Wind component m the meridional direction : 
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Temperature equation: 
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Pressure tendency equation: 

3P* f i 3(p*u) 3(p*v cos<j>) 

3t J o ^a cos <|>3A + "a cos 4> 3 ^ 

Moisture equation: 

3 (P*r) 3(p*ru) 3 (p*cos <J>rv) 3(op*r) 

3t a cos <f>3A a cos 4>3<j> 3a 


- C + 


F 

Hr 


+ 


F 

v r 


Hydrostatic equation: 


3$ _ RT 

Ikf c 


Vertical velocity (p) equation: 

3P* 3P* 9 p * 

w = p 0 + 0 [-— — + u rin — + v — vx] 

3t a cos <p3A a3<j> 

In the above equations , 
u = zonal wind component 
v = meridional wind component 
T = temperature 
r = moisture mixing ratio 
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p* = surface pressure 

$ = gz = geopotential 

f - Coriolis parameter 
* (3.cx 

o — = vertical velocity (a) 

ui = vertical velocity (p) 

Q = diabatic heating term 

C = moisture source/s ink term 

Cp ~ specific heat of air at constant pressure 

F = frictional force terms due to subgrid scale mixing 
in the horizontal and vertical directions 

The subscript p attached to the pressure gradient terms 

indicates differentiation on an isobaric surface. 

4 . 0 Numerics 

Finite difference methods are used to approximate the 
solution of the atmospheric model equations. 

4.1 Spacial Differencing 

The spacial finite differencing is based on the concept 
of the "box method" proposed by Bryan and is essentially 
similar to the energy conserving system suggested earlier 
by Lilly. This method basically consists of differencing 
quantities that have been averaged over one grid interval 
and is second-order . 
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The pressure gradient force terms m the momentum 
equations are evaluated on locally synthesized pressure 
surfaces to reduce the truncation error that arises from 
differencing these terms on sigma surfaces over terrain 
of varying height. 

Horizontal mixing of momentum, heat and moisture is 
based on the deformation scheme of smagonnsky with the 
characteristic scale of subgrid eddies defined as the 
mean of the meridional and Fourier filtered latitudinal 
grid distance. 

4.2 Temporal Differencing 

Centered time differencing is used with a time step 
of ten minutes for the N19 version of the model. 

The overspecification of tendencies due to the poleward 
convergence of the meridians is removed through Fourier 
filtering of all of the prognostic variables at each time 
step. The meridional component of the divergence operator 
is also Fourier filtered. 

The minimum wavelength allowed is specified by 


L MIN 


2 ira cos 4> c 



[D.9J 


where <f> is the equatorward boundary of the filtered region 
c 

and N is the number of grid points on a latitude circle. 
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(4> = 45° and N = 64 for the N19 version of the model.) 

c 

This minimum wavelength is obtained by limiting the 

maximum wave number at latitude > <f> to 

c 


_ 2Tra cos <f> _ N cos <J> 

*MAX L.. t „ 2 cos <P 

MIN c 


[D.103 


In the case of the vector variables u and v, the components 
are transformed to a polar stereographic projection before 
the Fourier filtering is performed. 

4 . 3 Boundary Conditions 

The horizontal boundary condition is periodicity in 
longitude and the poles are not grid points. The vertical 
boundary conditions are that o = 0 at a = 1 and 0 . 

5.0 Processes and Effects 

5.1 Dynamics 

Terms representing the horizontal and vertical advection 
of heat, moisture and momentum are included as are terms 
representing the pressure gradient and Coriolis forces 
m the momentum equations and a term representing adiabatic 
heating/cooling in the thermodynamic equation. 

Orography that is consistent with reality and the 
grid resolution is included to add realism to the model 
simulations . 
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Surface stress which is dependent on the surface wind 
and a terrain height dependent drag coefficient is 
computed. This surface stress constitutes the lower boundary 
condition for the computation of the Reynolds stress from 
vertical mixing. 

5 . 2 Heat Sources/Sinks 

Sensible heating m the lowest sigma layer is evaluated 
using a formulation for the heat flux at the surface of 
the earth that is similar to that used for the surface 
stress. The ground temperature over the continents is 
predicted from the net surface heating or cooling. The 
ocean temperature is spacially varying, but held constant 
during the simulation. 

Subgrid scale vertical mixing of heat is included in 
the uppermost eight layers. 

Solar and terrestrial radiation effects are included 
m the manner described by Manabe and Stnckler and Manabe 
and Wetherald and used in earlier GFDL models (Holloway, 
et al., 1971). The diurnal variation of the solar radiation 
is eliminated by the use of an effective mean zenith angle 
for each latitude. High, middle and low clouds are considered. 

Release of latent heat due to large scale precipitation 
is included. 

A moist convective adjustment procedure is included 
which has the net effects of stabilizing the lapse rate, 
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condensing water vapor, releasing the heat of condensation 
and transferring heat from lower to upper layers. 

5 . 3 Moisture Sources/Sinks 

The large scale precipitation algorithm removes 
supersaturated states from the model atmosphere which results 
m precipitation. 

The moist convection process can result m the 
condensation of moisture which also contributes to the 
precipitation . 

Evaporation from the surface (land or ocean) is evaluated 
using a formulation for the moisture flux that is similar 
to that used for the surface stress. 

Subgrid scale vertical mixing of moisture is included 
m the uppermost eight layers. 

6 . 0 Inputs 

The GFDL atmospheric general circulation model was 
developed for use m long range prediction and general 
circulation simulation experiments. Thus, data inputs 
are either self generated by going through a spin-up 
period for general circulation experiments or obtained 
from analyses of the atmosphere at a particular time for 
prediction experiments. These analyses, obtained from 
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various sources, are then interpolated to the model grid 
to obtain initial data for the forecast. 

7 . 0 Outputs 

An output cycle may be specified at any desired 
interval and produces the sea level pressure, heights and 
temperatures at standard levels, and numerous other 
forecast and diagnostic quantities. The fields on pressure 
surfaces are produced by interpolation from the sigma 
surface data. 


8 . 0 Computer Information 


The atmospheric general circulation model is programmed 
m Fortran IV for the Texas Instrument ASC Computer System. 


9 . 0 References 


Joint Organizing Committee - GARP, "Modelling for the First 
GARP Global Experiment," GARP Publication Series No. 

14, June 1974, Chapter 2, pp. 7-26. 

Holloway, J.L.Jr. and S. Manabe, "Simulation of Climate by 
a Global General Circulation Model," Monthly Weather 
Review, Vol . 99, No. 5, May 1971, pp. 335-370. 

Holloway, J.L.Jr., M.J. Spelman and S. Manabe, "Latitude- 

Longitude Grid Suitable for Numerical Time Integration 
of a Global Atmospheric Model," Monthly Weather Review, 
Vol. 101, No. 1, January 1973, pp. 69-78. 
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E. GFDL 18 -Layer Global Model 


The 18-layer Geophysical Fluid Dynamics Laboratory (GFDL) 
atmospheric general circulation model was developed for use 
m long range prediction and general circulation experiments . 
This model is a predecessor of the 9-layer model described 
in Section D and the governing equations and many of the 
parametenzations are identical. Thus, reference is made 
to Section D where appropriate rather than repeat the 
material here. Of particular interest in this model is 
the finite difference grid that is used - the modified 
Kurihara grid. 


MODEL DESCRIPTION 

1.0 Variables 

1 . 1 Independent 

The independent variables are the longitude, X, and 
the latitude, <j> ; sigma in the vertical defined as 



where p* is the surface pressure; and time, t. 
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1 . 2 Dependent 


1.2.1 Prognostic 

The prognostic variables are the longitudinal component, 
u, and the meridional component, v, of the wind; the 
temperature, T; the surface pressure, p*; and the moisture 
mixing ratio, q. Secondary prognostic variables are the 
surface temperature over land and the precipitation. 

1.2.2 Diagnostic 

The diagnostic variables are the geopotential , 4> , 
the vertical a velocity, a = , and the vertical p 

velocity, u. 

2 . 0 Domain 

2.1 Horizontal 

The horizontal grid is a global latitude-longitude 
grid with \ and <Js as coordinates . Grid points are placed 
on the globe according to a scheme known as the modified 
Kunhara grid. Figure E-l shows a portion of a modified 
Kurihara grid which has nineteen grid points between the 
equator and a pole. (There are 48 points between a pole 
and the equator m this model giving rise to the denotation 
of N48.) Variables are placed at the center of the grid 
boxes and, thus, the poles are not grid points. There is 
no staggering of variables m the horizontal. 
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Figure E-l 

Polar Stereographic Projection of a Portion of a 
Modif led Kurihara Grid 


2.2 Vertical 

The vertical coordinate is sigma as defined by 
Equation [E-l] . Figure E-2 shows the vertical structure 
of the model and the placement of the sigma surfaces 
and variables m the vertical. 
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Figure E-2 

Vertical Structure of the GFDL 18-Layer Global Model 


2.3 Temporal 

Two time levels of data are required (except during 
the startup step) since a centered time differencing scheme 
is used along with periodic Matsuno time steps . There is 
no staggering of variables m time. 
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3 . 0 Governing Equations 


Refer to Section D-3.0. 

4 . 0 Numerics 


Finite difference methods are used to approximate 
the solution of the atmospheric model equations. 

4.1 Spacial Differencing 

The spacial finite differencing is based on the concept 
of the "box method" proposed by Bryan. In the modified 
Kurihara grid each box is surrounded by six other boxes 
as shown m Figure E-3 and all variables are defined at 
the centers of the boxes. The horizontal flux 


N (y) 


3 (P*uy) 3 (p*vcos<f>y) 
acos<J>3A + acos$ 3 <j) 


[E.2] 


is approximated by the finite difference operator 


N (y) = ^ 


(P*V) + (P*V) 


+ y, 


w f ] 


i=i 


[E.3] 


where l is the index for the surrounding boxes, V is 
the horizontal flow vector, w is the box weight and the 
summation is made for the surrounding boxes. Other terms 
are approximated m a similar manner. 
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Figure E-3 

Modified Kunhara Grid Box Indexing Scheme 

The pressure gradient force terms in the momentum 
equations are evaluated on locally synthesized pressure 
surfaces to reduce the truncation error that arises 
from differencing these terms on sigma surfaces over 
terrain of varying height. 

Horizontal mixing of momentum, heat and moisture is 
based on the deformation formulation of nonlinear viscosity 
of Smagonnsky. 
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4.2 Temporal Differencing 


Centered time differencing is used along with periodic 
Matsuno time steps for control of computational noise . 

The overspecif ication of tendencies due to the 
poleward convergence of the meridians is removed either 
through a spectral filtering technique or weighted averaging 
at high latitudes (the six rows nearest each pole) . 

4 . 3 Boundary Conditions 

The horizontal boundary condition is periodicity m 
longitude and the poles are not grid points. The 
vertical boundary conditions are that cr = 0 at o = 1 and 0 . 

5.0 Processes and Effects 
Refer to Section D-5.0. 

6 . 0 Inputs 

The 18-layer GFDL atmospheric prediction model was 
developed for use in long range prediction and general 
circulation simulation experiments. Thus, data inputs 
are either self generated by going through a spin-up 
period for general circulation experiments or obtained 
from analyses of the atmosphere at a particular time 
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for prediction experiments. These analyses are then 
interpolated to the model grid to obtain initial data 
for the forecast. 


7 . 0 Outputs 

An output cycle may be specified at any desired 
interval and produces the sea level pressure, heights and 
temperatures at standard levels , and numerous other 
forecast and diagnostic quantities. The fields on pressure 
surfaces are produced by interpolation from the sigma 
surface data. 


8 . 0 Computer Information 

The atmospheric general circulation model is programmed 
in Fortran IV for the Texas Instruments ASC Computer System. 

9 . 0 References 


Joint Organizing Committee - GARP, "Modelling for the 
First GARP Global Experiment," GARP Publication 
Series No. 14, June 1974, Chapter 3, pp. 27-41. 

Miyakoda, K. , R.W. Moyer, H. Stambler, R.H. Clarke and 
R.F. Strickler, "A Prediction Experiment with a 
Global Model Made in the Kunhara Grid," Journal 
of Meteorological Society of Japan, Vol. 49, 1971, 
pp. 521-536. 
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F. GISS 9-Layer Global Model 


The GISS global atmospheric prediction model was 
developed as a medium-range forecast model for use m 
observing system simulation experiments, asynoptic data 
assimilation studies and experimental long-range forecasting 
(Sommerville, et al. , 1974). The GISS model shares the 
overall structure of the UCLA, three-level model (Arakawa, 
1972) and retains the sigma coordinate formulation of the 
governing equations, the Arakawa numerical method with 
advective quasi-conservation of quadratic quantities and 
much of the UCLA representation of physical processes 
occurring at and near the lower boundary of the atmosphere. 
The GISS model differs from the UCLA model m having a 
finer vertical resolution {nine levels) and in its treatment 
of four areas of physical processes: moist convection, 

turbulent subgrid-scale processes, solar radiation and 
terrestrial radiation. 

MODEL DESCRIPTION 

1.0 Variables 

1 . 1 Independent 

The independent variables are A and (f> of a spherical 
polar coordinate system in the horizontal; sigma m the 
vertical defined as 
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a 


[F.l] 



where P is the pressure at the top of the model atmosphere 
and P is the pressure at the bottom of the model atmosphere; 
and time, t. 

1 . 2 Dependent 

1.2.1 Prognostic 

The prognostic variables are the pressure, tt = P - P , 
the temperature, T, the latitudinal wind component, u, and 
the longitudinal wind component, v, and the specific humidity, 
q. A secondary prognostic variable is the precipitation 
obtained from large scale precipitation and the moist 
convection process. 

1.2.2 Diagnostic 

The diagnostic variables are the geopotential , 

* dcj 

and the vertical velocity, o = — . 

2 . 0 Domain 

2.1 Horizontal 

The horizontal domain is a global latitude, longitude 
grid with 1 and <f> as the coordinates and the grid resolution 
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is 4° in latitude and 5° m longitude. The poles are 
considered as mass points and are treated separately. 

The variables are distributed m the horizontal according 
to Arakawa's Scheme B as shown m Figure F-l. 
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Figure F-l 

Horizontal Grid Structure of the GISS 9-Level Global Model 


2.2 Vertical 

The vertical coordinate is sigma as defined by Equation 
[F . 1] . Figure F-2 shows the model sigma surfaces and the 
placement of the variables m the vertical. The vertical 
levels are equally spaced in sigma. 
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Figure F-2 

Vertical Structure of the GISS 9-Level Global Model 
2 . 3 Temporal 

One time level of data is required since either the 
Time Averaged Space Uncentered {TASU) - Matsuno or regular 
Matsuno time steps are used. Each simulated half-hour 
consists of ten minutes of TASU - Matsuno steps and 
twenty minutes of regular Matsuno steps with a time step 
of five minutes. There is no staggering of variables m 
time . 
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3 . 0 Governing Equations 

The governing partial differential equations written 
m flux form with a sigma coordinate m the vertical are 
listed below. 

Equation of motion: 


+ fkxV + V# + craViT = F 
Pressure tendency equation: 

H = -jf ’•(irviaa 

First law of thermodynamics : 

de _ Q 6 

dt C T 
P 

Water vapor equation: 


[F.2] 


[F.3] 


[F.4] 



-C + E 


Equation of state: 


[F.5] 


Pa = RT 


[F.61 



Hydrostatic equation: 


V 
f 
k 

V 
a 

Tf 

F 

R 

T 

e 

c 

] 

Q 

$ 

q 

c 

E 


= -air 
3a 


[F.7] 


In the above equations : 

= horizontal velocity 
= Coriolis parameter 
= vertical unit vector 
— two dimensional gradient operator 
= specific volume 



= horizontal frictional force 
= gas constant 
= temperature 
= potential temperature 

= specific heat of dry air at constant pressure 

= heating rate per unit mass 
= gz = geopotential 
= water vapor mixing ratio 
= rate of condensation 
= rate of evaporation 
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4.0 Numerics 


Finite difference methods are used to approximate 
the solution of the forecast equations. 

4.1 Spacial Differencing 

The horizontal differencing follows that of Arakawa 
(1972} . The distribution of variables on the horizontal 
grid corresponds to his Scheme B as shown m Figure F-l. 

The space differencing of the nonlinear advective terms m 
the equation of motion is constructed so as to maintain a 
constraint analogous to the conservation of mean square 
vorticity. The TASU-Matsuno time steps use uncentered 
differencing while the Matsuno time steps use centered 
differencing . 

The polar calculations follow Arakawa m which the 
change of the variables defined at the poles (tt, T and q) 
is a result of the meridional flux at all points on the 
latitude circle around the particular pole. 

4 . 2 Temporal Differencing 

Time differencing proceeds m half-hour cycles using 
a five minute time step. The first two steps of a cycle 
are TASU-Matsuno and the last four are regular Matsuno steps 
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(Arakawa, 1972) . For computational economy, physical 
processes other than advection are calculated every 
half hour of simulated time, except for solar and 
terrestrial radiation calculations which are performed 
every two hours. 

To avoid reducing the time step due to the convergence 
of the grid near the poles, the zonal mass flux and zonal 
pressure force are smoothed longitudinally using the 
three point smoothing technique described by Gates et al . 

(1971) . 

4 . 3 Boundary Conditions 

The horizontal boundary conditions are periodicity m 
longitude and the poles are treated as mass points as 
explained m Section 4.1. 

The vertical boundary conditions are a = 0 at a ~ 0 
and 1, and the upper pressure P^_ is taken to be 10 mb. 

5.0 Processes and Effects 

5 . 1 Dynamics 

Terms representing the horizontal and vertical advection 
and convergence/divergence of heat, moisture and momentum are 
included as are terms representing the pressure gradient and 
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Coriolis forces m the equation of motion. Adiabatic 
temperature changes are also modeled. 

Orography that is consistent with reality and the grid 
resolution is included to add realism to the flows. 

The surface stress is parameterized by drag laws with 
the drag coefficient being a function of wind speed, type 
of surface and surface height (to simulate mountain drag) . 

The surface velocity is estimated by extrapolating the 
velocities in the lowest two layers to the lower boundary. 

5.2 Heat Sources/Sinks 

Sensible heating in the boundary layer (lowest sigma 
layer) is evaluated by the method of Arakawa (1972) with 
the diffusivity given by Deardorff used to evaluate the 
thermal flux. The ground temperature (land, ice or snow) 
is predicted from the net surface heating or cooling. The 
ocean temperature is spacially varying, but held constant 
during the forecast. 

Subgrid scale vertical diffusion of heat with a constant 
diffusivity is included in the upper eight sigma layers. 

Solar and terrestrial radiation effects are included 
as fairly detailed parameterizations (Sommerville, et al., 
1974.) Cloudiness effects are considered at all sigma levels. 

Release of latent heat due to large scale precipitation 
is also included. 
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A parameterization of cumulus ensembles (moist 
convection) based on the Arakawa (1972) method redistributes 
heat and moisture in the vertical. The lower six model 
layers are "strapped" into three layers analogous to the 
UCLA 3-layer model. A dry convective adjustment procedure 
is included to preclude hydrostatic instability. 

5.3 Moisture Sources/Sinks 

The large scale precipitation algorithm removes 
supersaturated states from the model atmosphere and may result 
m evaporation m layers below the saturated layer. 

The moist convection process can result in the 
evaporation/condensation of moisture m the lower six 
layers . 

Evaporation/condensation from the surface into the 
lowest sigma layer is evaluated as part of the Arakawa 
type boundary layer computations. Subgrid scale vertical 
diffusion of moisture with a constant diffusivity is 
included m the upper eight sigma layers. 

6.0 Inputs 

The GISS model while generally used as a forecast 
model, has also been tested as a climate simulation model 
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{long term January) . Thus, data inputs are either self 
generated by going through a spin-up period for climate 
simulation or obtained from analyses of the atmosphere at 
a particular time. These analyses, obtained from various 
sources, are then interpolated to the model grid and 
no further balancing or smoothing procedures are 
employed. 

7 . 0 Outputs 

An output cycle may be specified at any desired 
interval and produces the sea level pressure, heights and 
temperatures at standard levels, and numerous other 
forecast and diagnostic quantities. The fields on pressure 
surfaces are produced by interpolation from the sigma 
surface data. 

8 . 0 Computer Information 

The forecast model is programmed in Fortran IV 
{Level H Compiler) for the IBM 360/95 Computer System. 

The code is quite modular and the execution time is 70 
minutes of CPU time per simulated model day. 

9 . 0 References 
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G. NCAR Multi-Layer Second-Generation Global Model 

Development of the NCAR second-generation model was 
begun m 1968 as a rewrite of the original six layer, 5° 
global model m order to obtain efficient programming with 
some degree of modularity. Various changes in the physical 
processes were also incorporated at the time. Over the years, 
this model has become a set of models with various resolutions 
and physical processes. The basic versions may be categorized 
as : 

• A six layer (tropospheric) , 5° model 

• A twelve layer (stratospheric) , 5° model 

• A six layer (tropospheric) , 2 1/2° high resolution 
model . 

Thus, by necessity, the description that follows is 
somewhat vague concerning the exact physical processes con- 
tained m a given model version since any particular version 
has mutated over the years . 


MODEL DESCRIPTION 

1.0 Variables 

1.1 Independent 

The independent variables are the longitude, the 
latitude, <J), the geometric height, Z, and time, t. 
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1.2 Dependent 


1.2.1 Prognostic 

The prognostic variables are the pressure, p, at the 

various levels, the longitudinal component, u, and the 

meridional component, v, of the wind, and the specific 

humidity, q. Secondary prognostic variables are the ground 

temperature, T , and the soil moisture m certain versions. 

9 

1.2.2 Diagnostic 

The main diagnostic variables are density, P, temperature, 
T, and the vertical velocity, w. 

2 . 0 Domain 

2.1 Horizontal 

The horizontal domain is a global latitude, longitude 
grid with X and $ as the coordinates. The resolution is 
either 5° or 2 1/2° giving a 72 x 37 grid or 144 x 73 grid, 
respectively. The variables are staggered in space and time 
as shown m Figure G-l. The poles are treated separately 
and are considered as mass points for pressure. 
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O denotes time level NAt 

X denotes time level {N+l)At 
or ( N-l )At 


Figure G-l 

Horizontal Staggering of Variables m Space and Time 
for the Second-Generation NCAR Model. 
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2.2 Vertical 


The vertical coordinate is geometric height, Z. The 
vertical grid increment is 3 Km giving a model atmosphere top 
of 18 Km for the six layer (tropospheric) version and 36 Km 
for the twelve layer (stratospheric) version. Figure G-2 
shows the vertical grid and the placement of variables m the 
vertical for the six layer tropospheric version. 
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Vertical Grid and Placement of Variables for the 
Tropospheric Second-Generation NCAR Model. 
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2 . 3 Temporal 


Due to the staggering of the variables m space and time 
as shown in Figure G-l, only one copy of the history variables 
is required for the centered time difference scheme. A time 
step of six minutes is used with the 5° versions and a time 
step of three minutes is used with the 2 1/2° version. 


3 . 0 Governing Equations 

The governing partial differential equations written 
m spherical polar coordinates in the horizontal and geo- 
metric height in the vertical are listed below. 

Momentum equation in the longitudinal direction: 


3 (pu) _ _ 1 . 5 (puu) , 9 (uvpcos<{>) , _ 9 (pwu) 

3 1 acos<j) L 9A 9<}> J 92 


+ ( f + Hiasi, pv . 

a 


_J= l£ + F 

acos<{> 9 A A 


[G.l] 


Momentum equation m the latitudinal direction: 


9 (pv) _ _ 1 r 3 (puv) 9 (wpcoscj)) _ 9 (pwv) " 

3t ~ acos^i 1 3A 9<j> J 9Z 


(f + pu - I |E , , 


[G.2] 
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Moisture Conservation equation: 


3 (pq) _ _ 1 (pug) 3 (pvqcostj)) , _ 3 (pwq) 

9t acos<j) 1 3A 94> J az 


+ M + pE 


[G.3] 


Pressure tendency equation: 


IP 

at 


B + gpw 


-J; 


1 r 3 jpu) , 9 (pvcos4)) 1 
acos^ L 8A 34> J 


[G.4] 
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C p T acos4> l aA 34> J i 


- f Z T dZ 

Y 4 p 
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where 


j = u 3p . v 3p „ / Z T 

acos4> 8 A a 34) g "'z 


1 r 3 (pu) 3 (pvcos4) ) , 

acos4> L 3A 34> J 

[G. 6] 


and w = V • VH at Z = H. 
s s 


Density is evaluated using the hydrostatic equation: 


_ _ 
3Z 


pg 


[G.7] 


Temperature is computed using 


T = 


R gg (£np) 


[G . 8 ] 
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Vertical velocity is evaluated using Richardson's 


equation: 



1 

acos<J> 


r 3u , 3 (vcoscj>) , 
l 9A 3<j> J 
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In the above 
u = 

v = 

p 

q 

P 

p w = 

H 

T 


w 
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g 
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equations , 

zonal wind component 
meridional wind component 
pressure 

P w /p = specific humidity 

density of air 

density of water vapor 

terrain height 

temperature 

vertical velocity 

rate of heating per unit mass 

longitudinal and latitudinal components 
of frictional force per unit volume 

condensation rate of water vapor per unit 
volume 

rate of change of water vapor per unit 
mass due to the vertical and horizontal 
diffusion of water vapor 

radius of the earth 

acceleration due to gravity 

Coriolis parameter 

gas constant for dry air 

specific heat for dry air at constant 
pressure 

specific heat for dry air at constant volume 



4 . 0 Numerics 


Finite difference methods are used to approximate the 
solution of the primitive equations. 

4.1 Spacial Differencing 

Second-order centered differencing is used to calculate 
the horizontal derivatives and vertical derivatives except 
near the surface where first-order approximations are used. 

Grid points are dropped near the poles to maintain a 
nearly constant geographical distance between mesh points 
for linear stability. Table G-l defines the number of grid 
points per latitude circle of the 5° model. The points are 
placed on the latitude circles such that there is a point at 
X = - ir (180 °W) and the points on each line are equidistant 
(see Figure G-3) . When neighbors on the latitude circles 
A(f> above and below a point are needed they are linearly 
interpolated from the nearest two points on the appropriate 
circle at the appropriate level in time. This procedure, of 
course, renders the finite differencing first-order near the 
poles . 

At grid points neighboring terrain which penetrates 

some AZ level, the second-order differencing scheme must 

\ 

be modified somehwat. Rather than use special difference 
formulas at these points, virtual points are created by 
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extrapolation through the midpoint of the difference sten- 
cil. This allows the use of the second-order difference 
formulas, but since the virtual points were created by 
linear extrapolation, the computation becomes first-order 
near the blocking terrain. 

The poles are treated as mass points for pressure m the 
computation with the polar pressures being defined as the 
average of the pressures on the latitude circle surrounding 
the pole. 

The horizontal eddy viscosity terms are based on either 
the deformation formulation of Smagormsky or the divergence 
of the vorticity as formulated by Leith according to the 
particular model version. In both cases time levels are 
mixed, thereby avoiding grid separation which is a charac- 
teristic of centered time differencing. 



Figure G-3 


Grid Near The Poles for the NCAR Second-Generation 
Model . 


G-9 



Table G-l 


Number of Grid Points as a Function of Latitude 
for NCAR Second-Generation Model. 
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4.2 Temporal Differencing 


Centered time differencing is used with either a six 
or three minute time step for the 5° and 2 1/2° models, 
respectively . 

A non-linear pressure smoother is applied on every time 
step which mixes the time levels for additional computational 
stability. 

A characteristic of the second-generation model which 
is due to the interpolation schemes used near the poles and 
around domains blocked by the orography, is an increase in 
the globally averaged value of surface pressure with time. 
Therefore, a procedure is invoked to maintain the globally 
averaged value of surface pressure at each time step. 

4 . 3 Boundary conditions 

The horizontal boundary conditions are periodicity in 

longitude and the poles are treated as mass points for 

pressure as explained m Section 4.1. 

The vertical boundary conditions are w = 0 at Z = 

and w=w = V ’VH at Z = H. 
s s 

5.0 Processes and Effects 

5 . 1 Dynamics 

Terms representing the horizontal and vertical advection 
and convergence/divergence of momentum and moisture are 
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included as are terms representing the pressure gradient 
and Coriolis forces in the momentum equations . The hori- 
zontal and vertical diffusion of water vapor is included 
in the moisture equation, and eddy diffusion of heat is 
included m the pressure equation. 

Orography that is consistent with reality and the grid 
resolution of the particular version of the model is included 
to add realism to the flows. 

The frictional force includes horizontal and vertical 
(Deardorff) diffusion of momentum. Various formulations 
have been used for the horizontal diffusion including the 
deformation of Smagormsky and Leith's divergence of vorticity 
procedures . 

5.2 Heat Sources/Sinks 

Both solar and terrestrial radiation are included accord- 
ing to the formulation of Sasamori. Clouds are included 
at two levels (3 Km and 9 Km) m the tropospheric version. 

The stratospheric version includes heating due to ozone 
absorption of solar energy and numerous levels of cloudiness. 

The release of latent heat due to the condensation of 
water vapor is included as part of the precipitation algorithm. 

Heating/cooling due to eddy diffusion m both the hori- 
zontal and vertical (Deardorff) is included. 

A surface energy balance equation is used to determine 
T over land areas. 

g 
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A dry convective adjustment procedure is included to 
preclude hydrostatic instability. 

The high resolution model (2 1/2°) includes elaborate 
parameterization schemes which include a convective para- 
meterization (Krishnamurti and Moxim) . 

5.3 Moisture Sources/Sinks 

The precipitation algorithm removes supersaturated 
states with upward vertical motion from the model atmosphere 
as precipitation. 

Evaporation or condensation is included in the Ekman 
layer (Deardorff) as is the horizontal diffusion of water 
vapor . 

The high resolution model (2 1/2°) includes elaborate 

parameterization schemes for convective parameterization 

(Khrishnamurti and Moxim) and a predictive equation for snow 

* 

cover and soil moisture. 

6 . 0 Inputs 

The models are run m two basic modes : 

a) The "GCM mode" m which data is self generated 
by going through a spin-up period which starts with an 
isothermal atmosphere at rest to obtain climate simulations, 

b) the "real data mode" m which analyzed data sets 
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(either subjective or objective analyses) obtained from 
various sources are used to make short and mid-range fore- 
casts. 

7 . 0 Outputs 

An output cycle for the models consists of dumping the 
current integration state onto a history tape at regular 
intervals of model time. A large, sophisticated processor 
program then processes these history tapes to produce a 
wide variety of output fields which are under the control 
of the user. Numerous other processor programs which pro- 
duce quantities such as climate statistics, etc. exist 
to process the history tapes . 

8.0 Computer Information 

The models are programmed in NCAR Fortran (a dialect 
of FTN IV) for the CDC 7600 using the NCAR Operating System. 
Running times m CPU hr/day for the various versions are 
listed below: 


6 

layer. 

5° 

0.2 

hr. 

12 

layer. 

5° 

0.4 

hr. 

6 

layer , 

2.5° 

1.6 

hr. 
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H. NCAR Multi-Layer Limited Area Model 


The NCAR limited-area model (LAM) was developed for 
numerical weather prediction studies on a very fine scale. 
This model is a version of the second-generation model des- 
cribed in Section G and contains physical parametenzations 
identical to the low resolution (5°) model. Thus, m the 
description that follows , reference is made to Section G 
where appropriate rather than repeating the material m 
this section. 

MODEL DESCRIPTION 

1.0 Variables 

1.1 Independent 

The independent variables are the longitude, the 
latitude, <j>, the geometric height, Z, and the time, t. 

1.2 Dependent 
1.2.1 Prognostic 

The prognostic variables are the pressure, p, at 
the various levels, the longitudinal component, u, and the 
meridional component, v, of the wind, and the specific 
humidity, q. A secondary prognostic variable is the ground 
temperature , . 
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1.2.2 Diagnostic 


The main diagnostic variables are density, p , tempera- 
ture, T, and the vertical velocity, w. 

2 . 0 Domain 

2.1 Horizontal 

The horizontal domain is a limited area of the globe 
which depends on the horizontal resolution chosen (5° to 
5/8° grid increments) . Typical fine scale (5/8° grid incre- 
ment) forecasts are made for regions encompassing 90° of 
longitude and 55° of latitude. Figure G-l m Section G shows 
the placement of variables in the horizontal. 

2.2 Vertical 

The vertical coordinate is geometric height, Z. The 
vertical resolution is variable from six to twenty-four 
layers with various vertical grid increments. Figure G-2 
in Section G shows the placement of the variables m the 
vertical. 

2 . 3 Temporal 

The staggering of variables m space and time is shown 
m Figure G-l of Section G. The time step used depends on 
the horizontal resolution, ranging from six minutes for 5° 
resolution to three-quarters of a minute for 5/8° resolution. 
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3 . 0 Governing Equations 


Refer to Section G-3.0. 

4 . 0 Numerics 

Finite difference methods are used to approximate the 
solution of the primitive equations. 

4.1 Spacial Differencing 

A discussion of the spacial differencing techniques and 
the dropping of grid points near the poles is given in 
Section G-4.1. The areas near the poles (north or south 
of 60°) are not allowed m the LAM domain due to the problem 
arising from the linear interpolation required to obtain 
inflow boundary values from the global forecast. 

Increased diffusion is used near the edges of the LAM 
to suppress computational noise arising from the boundary 
calculations. The diffusion coefficient is increased by a 
factor of four m the 5° band along the boundary and by a 
factor of two m the adjacent 5° band. 

4.2 Temporal Differencing 

A discussion of the temporal differencing is given m 
Section G-4.2. No correction to the surface pressure is made 
m the LAM during the integration. 
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4 . 3 Boundary Conditions 


The vertical boundary conditions are w = 0 at Z = and 

w = w = V *VH at Z=H. 
s s 

The horizontal boundary values are determined from 5° 
global forecast values for the inflow boundary points and by 
Lagrangian extrapolation of the LAM interior point values 
for the outflow boundary points. The decision as to which 
boundary points are inflow points is based on the winds 
obtained from the 5° global forecast. If the winds are 
directed into the LAM domain, then the boundary point is an 
inflow point and the LAM prognostic variables are determined 
from the corresponding global values by linear space- 
time interpolation. If the global winds on the LAM boundary 
are directed out of the LAM domain, then the point is an 
outflow point and the LAM prognostic variables are deter- 
mined by an extrapolation from the LAM interior to the boundary. 
The extrapolation is simply the horizontal advection of LAM 
variables by the global wind. 

5.0 Processes and Effects 
Refer to Section G-5.0. 

6.0 Inputs 

Input to the LAM consists of an initial state obtained 
from either a general circulation model simulation computation 
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or through an analysis of real data. Also required are 
computed values for the boundary points generated by the 
global model for the "forecast" period. The "global" winds 
are used to determine whether the boundary points are inflow 
or outflow points. The inflow points are then determined 
from the global values by linear space-time interpolation and 
the outflow points by Lagrangian extrapolation from the LAM 
interior. 

\ 

7 . 0 Outputs 

An output cycle for the model consists of dumping the 
current integration cycle onto a history tape at regular 
intervals of model time. A large , sophisticated processor 
program then processes these history tapes to produce a wide 
variety of output fields which are under the control of the 
user. 

8 . 0 Computer Information 

The LAM is programmed m NCAR Fortran (a dialect of FTN IV) 
for the CDC 7600 using the NCAR Operating System. Running 
times in CPU hr/day for various versions are listed below 
(ranges are given since running time depends on the size of 
the forecast domain) . 
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6 

layer , 

1.25° 

1.0 

- 2.0 

hr 

12 

layer , 

1.25° 

2.0 

- 4.0 

hr 

24 

layer , 

1.25° 

4.0 

- 8.0 

hr 

24 

layer , 

0.625° 

32.0 

-64.0 

hr 
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I. FfCAR Multi-Layer Third-Generation Global Model 

Development of the NCAR third-generation model was 
begun m 1972 to obtain an improved and more flexible model. 
The procedure of dropping points near the poles and of 
mountain blocking of layers m the second-generation model 
resulted m first-order difference approximations m these 
regions. Williamson and Browning (1973) showed that substan- 
tial improvement could be made by restoring the skipped 
gndpomts and removing computationally unstable modes by 
Fourier filtering. Fourth-order differencing is also in- 
cluded to reduce truncation errors. The model uses a new 

vertical coordinate, transformed height, developed by Kasahara 

/ 

(1974) . Both the horizontal and vertical resolution is varia- 
ble. 

The physical processes are essentially the same as m 
the second-generation model. However, solar and infrared 
radiation calculations are improved and cloudiness can be 
included at all levels. A new boundary layer treatment is 
be mg deve loped . 
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MODEL DESCRIPTION 


1.0 Variables 

1 . 1 Independent 

The independent variables are the longitude, X, the la- 
titude, <j> ; the transformed geometric height, Z, given by 



where Z T is the height of the model top and H is the 
terrain height; and time, t. 

1 . 2 Dependent 

1.2.1 Prognostic 

The prognostic variables are the pressure, p, at the 
various levels, the longitudinal component, u, and the 
meridional component, v, of the wind, and the specific 
humidity, q. A secondary prognostic variable is the ground 
temperature , T g . ’ 

1.2.2 Diagnostic 

The mam diagnostic variables are density, p, temperature, 

dZ 

T, and the vertical velocity, w = ^ . 
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2 . 0 Domain 


2.1 Horizontal 

The horizontal domain is a global latitude, longitude 
grid with <j> and X as the coordinates. The resolution is 
either 5.6° or 2.8° giving a 64 x 33 or 128 x 65 grid, 
respectively. The poles are treated separately to maintain 
the fourth-order accuracy of the finite difference approxi- 
mations . There is no staggering of variables in the hori- 
zontal . 

2 . 2 Vertical 

The vertical coordinate is transformed geometric height, 
Z. Figure 1-1 shows the vertical grid and the placement of 
variables in the vertical. 

2 . 3 Temporal 

A centered difference scheme is used with a time step 
of either five or two and one-half for the 5.6° and 2.8° 
grids, respectively. Fourier filtering of the prognostic 
variables at each time step is performed to remove the 
computationally unstable modes caused by the convergence 
of the grid near the poles . 




i“l fD i/D 1+1/3 

► 

X 


Figure 1-1 

Vertical Grid and Placement of Variables for the 
Third-Generation NCAR Model. 


3 . 0 Governing Equations 

The governing partial differential equations written 
m spherical polar coordinates m the horizontal and trans- 
formed height m the vertical are listed below. 

Longitudinal equation 

3u „ 1 r 9 (uu) . 

3t acosij) L 3X 

+ (f + 2*22*)v - 
a. 


of motion: 


3 (uvcos<j>) . _ 3 (wu) 


3 4> ' 


32 


1 r 3$ , RT 3p, . _ 

acos<J> p 3A F x 


[1.2] 


Latitudinal equation of motion: 


3v 

3t 


acos^ 3X 

3<j) 

(f + utan *)u - 
a 

1 r 3$ 
a L 3cf> 


-] - 


3 (wv) 
3Z 


P 3<j>’ 

Moisture conservation equation: 


ia 

3 1 


1 r 3 (uq) , 3(vqcos<p)n 
TrrTT L- yr _ + J ~ 


acos<j> SX 


3 4> 


3 (wq) 
32 


[1.3] 


+ ^ M + E 


[1.4] 
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Pressure tendency equation: 


|| = B - gpw (Z T — H) 




B = |£. at Z = 0, and it is assumed that w = 0 at Z 
at 


- 0 , 


and , thus , 


w s + ( z 


t _ h) X 


B = 


YP 


Q , 1 r 3u , 3 (vcos(j> ) ii jg 

CpT acostj) L 3A 3<j> 


Z t ~H r 

~T 


r°d| 

Ji p 


[ 1 . 6 ] 


where 


J = 


u 3p v 3p 


acos<}> 3A a 3<j> 


+ g(Zrn-H) 


f! 1 rliMl . 3 (pucosfr) ldS ri 71 
T J Z acos<fi — 3A 3(f) 


and w = - V • V ( Z -H) at Z = 0. 

S S 1 


Density is evaluated using the hydrostatic equation: 


— = gp (Z -H) 

3 Z 


[ 1 . 8 ] 


Temperature is computed using: 


T = 


R ( Z -H) - — (£np) 
3Z 


;i.9] 
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Vertical velocity is evaluated using Richardson's equation: 


w 


dz (\ 1 r 3u . 3(vcos<j>)-, , B+J 

at = w s ' K { i^i* + — 1 + yp - 


Q 


C T 
P 


> dZ 


[I. 10] 


In the above 
u = 

v = 

p 

g 

P 
P 
H 
T 
W 

Q 


w 




M 


a 

g 

f 

R 

$ 


equations : 

zonal wind component 
meridional wind component 
pressure 

p / p = specific humidity 

density of air 

density of water vapor 

terrain height 

temperature 

dZ 

= vertical-velocity 

rate of heating per unit mass 

longitudinal and latitudinal components 
of frictional force per unit volume 

condensation rate of water vapor per 
unit volume 

rate of change of water vapor per unit 
mass due to the vertical and horizontal 
diffusion of water vapor 

radius of the earth 

acceleration due to gravity 

Coriolis parameter 

gas constant for dry air 

gz = geopotential 
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= specific heat for dry air at constant 
pressure 

= specific heat for dry air at constant 
volume 


■ V c 


v 


4 . 0 Numerics 


Finite difference methods are used to approximate the 
solution of the primitive equations. 

4.1 Spacial Differencing 

Fourth-order centered differences are used to calculate 
horizontal derivatives and second-order centered differences 
are used to calculate vertical derivatives. An option exists 
m the model to use second-order differencing in the hori- 
zontal. 

The poles are handled separately to maintain the fourth- 
order difference approximations by transforming to polar 
stereographic coordinates to perform the pole calculations. * 

Various horizontal diffusion terms are available. These 
include a spectral chopping technique m both the latitudinal 
and longitudinal directions m which all waves smaller than 
4 Ax are set to zero at a specified time interval, and a 
combination method which consists of fourth“order smoothing 
m latitude and spectral chopping m longitude. 
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4.2 Temporal Differencing 


Centered time differencing is used with either a five 
or two and one-half minute time step for the 5.6° and 2.8° 
models, respectively. 

The elimination of computationally unstable modes near 
the poles due to the convergence of the grid is accomplished 
by Fourier filtering the prognostic variables every time step. 

4 . 3 Boundary Conditions 

The horizontal boundary conditions are periodicity m 
longitude and the poles are handled separately to remove 
the coordinate system caused discontinuity and maintain 
fourth-order accuracy. 

The vertical boundary conditions are w = 0 at Z = 0 
and 1. 

5.0 Processes and Effects 

5 . 1 Dynamics 

Terms representing the horizontal and vertical advection 
and convergence/divergence of momentum and moisture are in- 
cluded as are terms representing the pressure gradient and 
Coriolis forces in the momentum operations. The horizontal 
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and vertical diffusion of water vapor is included in the 
moisture equation, and eddy diffusion of heat is included 
in the pressure equation. 

Orography that is consistent with reality and the grid 
resolution of the particular version of the model is included 
to add realism to the flows. Frictional forces are included. 

5.2 Heat Sources/Sinks 

Both solar and terrestrial radiation are included ac- 
cording to the formulation of Sasamon. Clouds may be 
included at all levels, but diagnostic cloudiness at two 
levels is currently used. 

The release of latent heat due to the condensation 
of water vapor is included as part of the precipitation 
algorithm. 

Heating/cooling due to eddy diffusion m both the hori- 
zontal and vertical (Deardorff) is included. 

A surface energy balance equation is used to determine 
T over land areas. 

g 

A dry convective adjustment procedure is included to 
preclude hydrostatic instability. 

5.3 Moisture Sources/Sinks 

i 

The precipitation algorithm removes supersaturated 
states with upward vertical motion from the model atmosphere 
as precipitation. 
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Evaporation or condensation is included m the boundary 
layer through the bulk transfer coefficient method (Deardorff ) . 
Horizontal diffusion of water vapor is also included. The 
ground is assumed saturated. 

6 . 0 Inputs 

The model is run m two basic modes; 

a) The "GCM mode" in which data is self generated by 
going through a spin-up period which starts with an isother- 
mal atmosphere at rest to obtain climate simulations. 

b) The "real data mode" m which analyzed data sets 
{either subjective or objective analyses) obtained from 
various sources are used to make short and mid-range forecasts . 

7. 0 Outputs 

An output cycle for the model consists of dumping the 
current integration state onto a history tape at regular 
intervals of model time. A large, sophisticated processor 
program then processes these history tapes to produce a wide 
variety of output fields which are under the control of the 
user. Numerous other processor programs which produce 
quantities such as climate statistics, etc. exist to process 
the history tapes. 

► 
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8.0 Computer Information 


The model is programmed xn NCAR Fortran (a dialect 
of FTN IV) for the CDC 7600 usxng the NCAR Operating 
System. Running times m CPU hr/day for various versions 


are listed below: 



6 layer. 

0 
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J. NMC 6 -Layer Hemispheric Operational Model 

The NMC operational atmospheric prediction model 
is currently executed twice daily to produce three day 
forecasts for the Northern Hemisphere. The model as 
described by Shuman and Hovermale (1968) was first put 
into operational use on 6 June 1966. Since that time 
numerous changes have been made to the model to improve 
the quality of the forecasts produced and to achieve a 
more efficient computation. The description that follows 
is, thus, a snapshot of the state of the prediction model 
as of mid 1976. 

MODEL DESCRIPTION 

1 . 0 Variables 

1 . 1 Independent 

The independent variables are x and y derived from a 
polar stereographic map projection m the horizontal; 
sigma in the vertical defined as : 



where P^ is the pressure at a given quasi-horizontal surface 
above some point and P T is the pressure at a surface below; 

Li 

and time , t . 
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1 . 2 Dependent 


1.2.1 Prognostic 

The prognostic variables are the x wind component, u; 

the y wind component, v; the potential temperature, 6; the 

vertical pressure gradient, ~ , at the three interfaces 

do 

(between the planetary boundary layer and the troposphere , 
between the troposphere and the stratosphere, and between 
the stratosphere and the constant 0 layer) ; and the specific 
humidity, q> m the planetary boundary layer and the tropo- 
sphere. A secondary prognostic variable is the precipitable 
water, W. 

1.2.2 Diagnostic 

The diagnostic variables are the geopotential, 

• da 

$, and the vertical velocity, a = ^ . 

2 . 0 Domain 

2.1 Horizontal 

The horizontal domain is a polar stereographic pro- 
jection (true at 60 °N) of the Northern Hemisphere using 
x-y coordinates with map factor 
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m 


1 + sm 60° 
1 + sm 4> 


[J.2] 


The grid resolution is 381 km at 60 °N and the grid 
dimensions are 53 x 57 which covers virtually the entire 
Northern Hemisphere as shown in Figure J-l. (The standard 
octagon is also shown.) There is no staggering of variables 
in the horizontal.- 

2.2. Vertical 

The vertical coordinate is sigma as defined by 
Equation [J-l] . Figure J-2 shows the vertical structure 
of the model with the four distinct sigma domains. The 
variables u, v, 0 and q are carried at the midpoints of 
the various sigma layers and the variables tt and p 
are carried at the layer interfaces. The vertical velocity, 
a T is defined at layer interfaces , but at the center of 
a horizontal gridpoint square. 
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Figure J-l 

NMC Northern Hemisphere Polar Stereographic Grid 
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Vertical Structure of the NMC Atmospheric Prediction Model 


2 . 3 Temporal 

Two time levels of data are required (except during 
the startup step) since a centered time differencing 
scheme is used with a twenty minute time step. There is 
no staggering of the variables m time. 
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3 . 0 Governing Equations 


The governing partial differential equations written 
for an x-y map projection with map factor m and vertical 
coordinate a are listed below. 

Wind component in the x direction : 


9u 
3 1 


- a 


3u 

3c 


+ v [ f - v ” + u — 3 + F + D 
ax ay x u 


[J.3] 


Wind component in the y direction: 


3v r 3v , av t 

TT-r- = -m [ u ~ — + V 2 — ] 
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3y p 3y 
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3 v 
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3x 3y y v 
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Thermodynamic energy equation: 


39 
3 1 


-m [ u + v &&• ] - o ~ + H 
3 x 3y 3c 
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Pressure tendency equation: 


3 t 
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Moisture equation: 


aa = 

St 


-m 




v 2a ] 
Sy J 


a aa + c 

Sa 


Hydrostatic equation: 


at. = -c 0 ax 

Sct ^p 0 c 


Exner function: 


* - ( raoo )K 

Xn the above equations , 

u = x direction wind component 

v = y direction wind component 

0 = potential temperature 

q = specific humidity 

$ = gz = geopotential 

p = pressure 

f - Coriolis parameter 

cr = ^2. = vertical velocity 

F x ,Fy = stress components 

D^,D v = lateral diffusion terms 

H = diabatic heatmg/cooling term 


[J.7] 


[J.8] 


[J.9] 
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C = condensation/evaporation term 

Cp = specific heat of dry axr at constant pressure 

R = gas constant for dry air 

< -b 

P 

4.0 Numerics 

Fxnxte dxfference methods are used to approximate 
the solutxon of the forecast model equatxons. 

4.1 Spacxal Dxfferencxng 

Centered, second-order fxnxte dxfferencxng xs used 
to approxxmate the horxzontal derxvatxves xn the various 
sxgma domaxns . A space averagxng technxque orxgxnated 
by Shuman (1968) is used in calculatxng the derxvatives 
which results in defining the derivatives at the midpoint 
of a grid square through averaging and then averaging 
these values back to the grid points. This averaging 
technique makes the finite difference scheme somewhat 
dissipative and, hence, more stable. 

The lateral diffusion terms included in the wind 
equations are of the form: 

D u = V-KVu [J.10] 
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where K = 0 for <j> > 20°N, K = 10 6 m 2 /sec for <p < 10 °N 
and K = 5 x 10 5 [1 + sm{ {t|» 0 — ) ] r <!> 0 = 15° and &<P - 5°. 

The vertical derivatives are approximated by second- 
order finite differences. 

4.2 Temporal Differencing 

Centered time differencing with pressure gradient 
force time averaging (Shuman, 1971) is used to obtain a 
time step of twenty minutes . The computation is started 
with an Euler backward or Matsuno time step. 

A light time smoothing of the wind and moisture 
fields is also performed. 

As a further control of the tropical region, the 
tendencies of all variables are truncated south of 20 °N. 

North of 20°N, the full tendency is used and south of 
10 °N one-half of the tendency is used. These two zones 
are blended using a latitudmally dependent function such 
that three-quarters of the tendency is used at 15 °N. 

4 . 3 Boundary Conditions 

The horizontal boundary conditions consist of insulated, 
slippery walls placed on the grid as shown in Figure J-l . 
That is, no heat or momentum flow is permitted through the 
wall, but a parallel flow is permitted along the wall. 



The vertical boundary condition that a = 0 for a 
material surface is specified at the top of the model, 
the boundary surface between the topmost computational 
layer and the stratosphere, the tropopause; and the 
earth's surface. As a consequence of the above approach, 
no mass exchange takes place among the lsentropic, 
stratospheric and tropospheric domains. However, mass 
exchange can take place between the boundary layer and 
the troposphere. 

5.0 Processes and Effects 

5 . 1 Dynamics 

Terms representing the horizontal and vertical 
advection and convergence/divergence of heat and momentum 
are included as are terms representing the pressure 
gradient and Coriolis forces in the momentum equations. 

Orography that is consistent with reality and the 
grid resolution is included to add realism to the flows. 

Horizontal stress is applied at the lowest level to 
simulate vertical momentum transport at the surface. 

5.2 Heat Sources/Sinks 

Both solar and terrestrial radiation are included. 
Since moisture is carried m the lower three layers only , 
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solar heating is evaluated m these layers. Cloudiness 
effects are considered in both the solar and terrestrial 
heating computations. 

Sensible heating is evaluated in the planetary boundary 
layer and a land-sea-ice discriminator is included. The 
ocean temperature is spacially varying, but held constant 
during the forecast. 

Release of latent heat due to large scale precipitation 
is included in the moisture bearing layers. 

A moist convective adjustment procedure is used to 
simulate sub-grid scale convective processes. Dry convec- 
tive adjustment to maintain hydrostatic stability occurs 
as a limiting case of the moist convective adjustment. 

5.3 Moisture Sources/Sinks 

Evaporation into the planetary boundary layer from 
the oceans is considered as a possible source of moisture. 

The moist convective adjustment procedure may 
redistribute moisture m the moisture carrying layers. 

The large scale precipitation algorithm acts as a 
moisture sink for precipitating layers. There is no 
absorption of precipitation in layers below precipitating 
layers . 
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6 . 0 Inputs 


The forecast model initial data is obtained from 
Cressman type hemispheric analyses of the mass structure 
(temperature and geopotentials) between 1000 mb and 100 mb, 
the sea level pressure, and the sea surface temperature. 
Analyses of the winds on the standard pressure surfaces 
are also made. 

The analyzed fields on pressure surfaces are then 
interpolated (linearly in £np) to obtain data on the 
model sigma surfaces. 

7 . 0 Outputs 

Temperatures , geopotential heights and winds on the 
mandatory pressure surfaces are produced by interpolation 
from the model sigma surfaces during an output cycle. 
Averaging of the output fields over time intervals 
surrounding the output time and filtering is used to 
render the output fields more visually pleasing. 
Precipitation fields are also produced during an output 
cycle. 

8 . 0 Computer Information 

The forecast model is run on the NOAA IBM 360/195 
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Computer System and is normally executed twice daily. 
The model is largely programmed m Fortran (Level H 
Compiler) with some specialized routines m assembly 
language and is fairly modular. 
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K. NMC 6 -Layer Limited Area Operational Model 

The NMC operational limited area, half-mesh model 
(LFM) is used m conjunction with the hemispheric model 
to obtain more detailed 48 hour forecasts over the North 
American Region. Boundary conditions for the LFM are 
obtained from the hemispheric model. The LFM is a 
modified version of the hemispheric model described m 
Section J. Thus, m the description that follows, 
reference is made to Section J where appropriate rather 
than repeating the material m this section. 


MODEL DESCRIPTION 

1 . 0 Variables 

1 . 1 Independent 

The independent variables are x and y derived from 
a polar stereographic map projection in the horizontal; 
sigma in the vertical defined as 



where Py is the pressure at a given quasi-horizontal 
surface above some point and P^ is the pressure at a 
surface below; and time, t. 
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1.2 Dependent 


1.2.1 Prognostic 

The prognostic variables are the x wind component, u; 
the y wind component, v; the potential temperature, 0; the 
vertical pressure gradient, , at the three interfaces 

<7 O 

(between the planetary boundary layer and the troposphere, 
between the troposphere and the stratosphere, and between 
the stratosphere and the constant 0 layer) ; and the specific 
humidity, q, m the planetary boundary layer and the tropo- 
sphere. A secondary prognostic variable is the precipitable 
water , W . 

1.2.2 Diagnostic 

The diagnostic variables are the geopotential, 

« cl(T 

and the vertical velocity, a = . 

2 . 0 Domain 

2.1 Horizontal 


The horizontal domain is a polar stereographic pro- 
jection (true at 60 °N) of part of the Northern Hemisphere 
using x-y coordinates with map factor 


1 -f sm 60° 
1 + sin <J) 


[K.2] 
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The grid resolution is 190.5 km at 6 0°N and the LFM grid 
covers the North American region. (See Figure J-l.) 

There is no staggering of variables m the horizontal. 

2 . 2 Vertical 

The vertical coordinate is sigma as defined by 
Equation [K.l] . Figure J-2 m Section J shows the vertical 
structure of the model with the four distinct sigma domains. 
The variables u, v, 0 and q are carried at the midpoints 
of the various sigma layers and the variables $, it and p 
are carried at the layer interfaces. The vertical velocity, 
a , is defined at layer interfaces , but at the center of a 
horizontal gndpomt square. 

2 . 3 Temporal 

Two time levels of data are required (except during 
the startup step) since a centered time differencing scheme 
is used with a six minute time step. There is no staggering 
of the variables in time. 

3 . 0 Governing Equations 
Refer to Section J-3.0. 

4 . 0 Numerics 

Finite difference methods are used to approximate the 
solution of the forecast model equations. 
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4.1 Spacial Differencing 


Centered, second-order finite differencing is used 
to approximate the horizontal derivatives in the various 
sigma domains. A space averaging technique originated 
by Shuman (1968) is used m calculating the derivatives 
which results m defining the derivatives at the midpoint 
of a grid square through averaging and then averaging 
these values back to the grid points. This averaging 
technique makes the finite difference scheme somewhat 
dissipative and, hence, more stable. 

The vertical derivatives are also approximated by 
second-order finite differences. 

4.2 Temporal Differencing 

Centered time differencing is used with a time step 
of six minutes . 

Two types of smoothers are used to control computational 
noise — spacial and temporal . A spacial smoother consisting 
of a twenty-five point operator is used at all points two 
or more grid intervals from the boundaries and a nine 
point operator is used at points one grid interval from 
the boundaries. The one-dimensional version of the 



twenty-five point operator is given by: 

^ = (1 - 6p )u x + 40(u i _ 1 + u 1+1 ) 

-B Z (Ui_ 2 + u x+2 ) [K.3] 

A value of 3 = 0.1 is currently used in the smoothing 
operator which is applied to all fields except precip- 
itable water. 

Robert time filtering is applied to all the fields — 

precipitable water, winds, and 0. 

3a 

4 . 3 Boundary Conditions 

The LFM horizontal boundary conditions are obtained 
from the hemispheric six-layer model six-hourly output 
fields. The tendencies on and near the boundaries consist 
of a blend of the LFM and coarse mesh tendencies given by: 

3X 

X n+l (I) = X n-l {1) + W{1) 3tf!l 2At 

9X J^ , 

+ [1 - w(l) 3 2 At ( [K.4] 

where X is the predicted variable, the subscript m denotes 
the LFM tendency and the subscript Li denotes the large- 
scale model tendency. The weighting factor W(I) is 
given by: 
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= 0 for I = boundary points 

= 0.33 for I = boundary points - 1 grid length 

W(I) 

= 0.67 for I = boundary points - 2 grid length 
= 1.00 for I = all other interior points 
The coarse mesh model tendencies are calculated from six- 
hourly history fields biquadratically interpolated m 
space to the LFM grid points. (See Perkey and Kreitzberg, 
1976.) 

The vertical boundary conditions are identical to 
those of the coarse mesh model as described m Section J-4 . 3 . 

5.0 Processes and Effects 

The various physical processes and effects are mostly 
described m Section J-5.0 since the LFM parametenzations 
differ only m detail from the six-layer hemispheric model . 

6 . 0 Inputs 

Initial data comes from a Cressman-type analysis on 
the LFM data points using raw data received during the one 
hour and forty-five minute time period after the synoptic 
time . 

The mass structure variables are initialized from 
geopotentials analyzed on the ten mandatory pressure levels; 
temperature at the earth's surface and at the mandatory 
pressure levels with the exception of 1000 mb; and the 
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pressure and temperature of the analyzed "tropopause" 
surface. These data fields are then interpolated to 
the sigma surfaces. 

Wind data is obtained from two sources, the coarse 
mesh twelve hour forecast winds on sigma surfaces', and the 
isobaric analysis of the winds interpolated to the sigma 
surfaces. A nondivergent wind is then calculated for each 
sigma layer which has the same relative vorticity as the 
interpolated analysis winds. Next, an irrotational wind is 
calculated which has the same divergence as the interpolated 
coarse-mesh forecast winds. The initial LFM wind is then 
taken to be the sum of the calculated nondivergent and 
irrotational wind components. 

7 . 0 Outputs 

Temperatures, geopotential heights and winds on the 
mandatory pressure surfaces are produced by interpolation 
and filtering from the model sigma surfaces during a six- 
hourly output cycle. Precipitation fields are also 

* 

produced during an output cycle. Various other parameters 
are also produced during an output cycle (e.g., MOS 
forecasts and input to aviation winds forecasts) . 

8 . 0 Computer Information 

The LFM forecast model is written m Fortran IV 
(Level H Compiler) , except for some specialized I/O 
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routines, for the NOAA IBM 360/195 Computer System. 

The code is fairly modular and requires 520 k hytes of 
memory and 25 million bytes of mass storage. The 
execution time is three and one-half minutes of CPU 
time per twelve hours of forecast. 
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L. NMC 8 -Layer Global Model 


Initial development of the NMC 8-layer atmospheric 
prediction model was completed in 1972. This model is an 
extension of the earlier 6-layer hemispheric model 
(described m Section J) having increased vertical 
resolution and a latitude- longitude grid allowing either 
a global or hemispheric domain. While the original 
version of this model embodied much of the physics of 
the older 6-layer model, many additions and improvements 
have been made which are outlined in the following 
description. 


MODEL DESCRIPTION 

1 . 0 Variables 

1 . 1 Independent 

The independent variables are the longitude. A, and 
the latitude, <j> ; sigma defined as 



where is the pressure at a given quasi-horizontal 
surface above some point and P L is the pressure at a 
surface below; and time, t. 
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1.2 Dependent 


1.2.1 Prognostic 

The prognostic variables are the longitudinal 
component, u, and the meridional component, v, of the 
wind; the potential temperature, 0; the vertical 
pressure gradient, at the three interfaces (between 
the planetary boundary layer and the troposphere, between 

i 

the troposphere and the stratosphere and between the strato 
sphere and the constant 0 layer) ; and the specific humidity 
q, in the planetary boundary layer and the troposphere. 
Secondary prognostic variables are the precipitable water, 
W, and the surface temperature, T^. 

1.2.2 Diagnostic 

The diagnostic variables are the geopotential , 

$, and the vertical velocity", o = ^ . 

2 . 0 Domain 

2.1 Horizontal 

The horizontal domain is a global (or hemispheric) 
latitude- longitude grid with X. and <p as the coordinates. 

The resolution is either 2 1/2° or 2° (with 5° as a test 
version) . The poles are treated separately as averages 
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of the forecast values one row removed from the pole . 
There xs no staggering of variables m the horizontal. 

2 . 2 Vertical 

The vertical coordinate is sigma as defined by 
Equation [L.l]. Figure L-l shows the vertical structure 
of the model with the four distinct sigma domains. The 
variables u, v, 6 and q are carried at the midpoints of 
the various sigma layers and the variables $, t f and p 
are carried at the layer interfaces. The vertical 
velocity, a, is defined at layer interfaces, but at the 
center of a horizontal gridpomt square. 

2 . 3 Temporal 

Two time levels of data are required {except during 
the startup step) since a centered time differencing 
scheme is used with a ten minute (2.5°) or seven and one- 
half (2°) minute time step. There is no staggering of 
variables m time. 
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Level: 


k 

+ 



-1 p 


0 


0 p * 50 mb 

1 

2 

3 p * 300 mb 

4 


5 


6 p =>= 903 mb 

7 

8 p *= 1000 mb 


Fxgure L-l 

Vertical Structure of the NMC 8-Layer Model 
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3 . 0 Governing Equations 

The governing partial differential equations written 
m spherical polar coordinates and a vertical coordinate 
sigma are listed below: 

Wind component m the latitudinal direction: 


3u 

at 



v 


IS - a §^ + v(f + U •) 


ay 


- - c e ^S-+F 

ax p ax x 


[L.2] 


Wind component m the meridional direction: 


- - u |Z . v |v _ • 3v 
3x 3y 3a 


- u (f + 


utanj) 


M. _ 

ay 


V 


In 

3 y 


+ F 


[L.3] 


Thermodynamic energy equation: 


30 = 36 

at 3x 



30. 

3a 


+ H 


IL.4] 
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Pressure tendency equation: 




= — — — u) 

3x ^ 3a 


3 

by 


‘I? v) 


3 

3a 


(IE a) 
3a 


, /3|\ v tan 4> 

3a a 


Moisture equation: 


ia=- u 22- v 3S--aa +E p 

at ax sy 3a 

Hydrostatic equation: 


M. = -c e 
3a P 3a 

Exner function: 


ir 


( 


1000 


') 


iC 


where 


3 

3X 


3 

a coS(j)3X 


and 



3 

a3<J> 


[L.5] 


[L.6] 
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[L.8] 


[L.9] 
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In the above equations , 


u 

v 

e 

q 

$ 

p 

f 

a 


8 A 

= a cos<j) -r-TT = zonal wind component 

<7 t 

= a = meridional wind component 


= potential temperature 
= specific humidity 
= gz = geopotential 
= pressure 

= Coriolis parameter 
= —• = vertical velocity 


F ,F 
x' y 


H 

EP 

C 


R 


stress components 

diabatic heating/cooling term 

evaporation and precipitation processes which alter q 
specific heat of dry air at constant pressure 

gas constant for dry air 


K 



4 . 0 Numerics 


Finite difference methods are used to approximate 
the solution of the forecast model equations. 


4 . 1 Spacial Differencing 

Centered, second-order finite differencing is used to 
approximate the horizontal derivatives m the various sigma 
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domains. A space averaging technique originated by- 
Shuman (1968) is used m calculating the derivatives 
which results m defining the derivatives at the midpoint 
of a grid square through averaging and then averaging 
these values back to the grid points. This averaging 
technique makes the finite difference technique somewhat 
dissipative and, hence, more stable. 

The vertical derivatives are approximated by second- 
order finite differences. 

The values of the variables at the poles are obtained 
by averaging the variables one row equatorward and assigning 
this average to the pole value. The one vector wind at 
the pole is resolved into the u and v components appropriate 
to the meridian along which the pole is approached. 

4.2, Temporal Differencing 

Centered time differencing is used with a time step of 
either ten or seven and one-half minutes for the 2.5° and 
2° versions, respectively. Robert time filtering with 
fairly heavy weighting (1/4, 1/2, 1/4) is used to maintain 
computational stability. The computation is started with a 
forward time step. 

The overspecif ication of tendencies due to the poleward 
convergence of the meridians is handled by an averaging 
procedure. The tendency terms are averaged using a triangular 



weighting function that expands longitudinally with 
increasing latitude. This weighting function is given by: 

w “ 1 ” L sec d> x = 1//2 ' 1 1 / 2 ' • • * L sec 4> 0 

T o 

=0 x ^ L sec <J) Q [L.10] 

where <j> Q is the latitude of the tendencies being averaged 
and L = 1.5. 

4 . 3 Boundary Conditions 

The horizontal boundary conditions consist of period- 
icity m longitude and the pole calculations as described 
above (Section 4.1) for the global versions of the model. 

The additional boundary condition of horizontal symmetry 
about the equator is used m the hemispheric versions . 

The vertical boundary condition that a = 0 for a 
material surface is specified at the top of the model, 
the boundary surface between the topmost computational 
layer and the stratosphere, the tropopause; and the 
earth's surface. As a consequence of the above approach, 
no mass exchange takes place among the lsentropic, 
stratospheric and tropospheric domains. However, mass 
exchange can take place between the boundary layer and 
the troposphere. 
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5.0 Processes and Effects 


5 . 1 Dynamics 

Terms representing the horizontal and vertical 
advection and convergence/divergence of heat and momentum 
are included as are terms representing the pressure 
gradient and Coriolis forces in the momentum equations . 

Orography that is consistent with reality and the 
grid resolution is included to add realism to the flows . 

Horizontal stress according to Shuman (1968) is 
applied at the lowest level to simulate vertical 
momentum transport at the surface. 

5 . 2 Heat Sources/Sinks 

Both solar and terrestrial radiation are included. 

Since moisture is carried in the lower five layers only, 
solar heating is evaluated m these layers. Cloudiness 
effects are considered in both the solar and terrestrial 
heating computations. 

Sensible heating is evaluated m the planetary boundary 
layer and a land-sea-ice discriminator is included. The 
ocean temperature is spacially varying, but held constant 
during the forecast while the land temperature is predicted 
through an energy balance equation. 
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Release of latent heat due to large scale precipitation 
is included in the moisture bearing layers. 

A moist convective adjustment procedure is used to 
simulate sub-grid scale convective processes. Dry convective 
adjustment to maintain hydrostatic stability occurs as 
a limiting case of the moist convective adjustment. 

5 . 3 Moisture Sources/Sinks 

Evaporation mta the planetary boundary layer either 
from the oceans or over land is considered as a possible 
source of moisture. 

The moist convective adjustment procedure may redis- 
tribute moisture in the moisture carrying layers. 

The large scale precipitation algorithm acts as a 
moisture sink for precipitating layers. There is absorp- 
tion of precipitation m layers below precipitating layers 
if they are not saturated. 

6 . 0 Inputs 

The forecast model initial data is obtained from Hough 
function analyses {Flattery, 1970) . The end results of 
these analyses are heights, winds and temperatures at the 
mandatory levels up to 50 mb plus surface temperature 
(land and sea) , tropopause pressure, snow/ice cover, and 
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relative humidity at mandatory levels up to 300 mb. 

The analyzed fields on pressure surfaces are 
then interpolated (linearly m -Cnp) to obtain initial 
data on the model sigma surfaces. 

7 . 0 Outputs 

Temperatures, geopotential heights and winds on the 
mandatory pressure surfaces are produced by interpolation 
from the model sigma surfaces during an output cycle. 
Various filters are used to render the fields more 
visually pleasing. Precipitation fields are also 
produced during an output cycle. 

8.0 Computer Information 

The forecast model is run on the NOAA IBM 360/185 
Computer System. The model is largely programmed m 
Fortran IV (Level H Compiler) except for assembly language 
coding of high usage sections and some specialized I/O 
routines. The code is fairly modular and requires 750 k 
bytes and 950 k bytes of storage for the 2.5° and 2° 
versions, respectively. The mass storage requirement is 
approximately 50 million bytes. The execution time is 
ten CPU minutes per twelve hour forecast for both the 2.5° 
global and 2° hemispheric models. 
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M. UCLA 12 -Layer Global Model 

The UCLA Mmtz-Arakawa atmospheric general circulation 
model was developed for numerical simulation experiments 
which study the mutual interactions of various physical 
processes m the atmosphere. This latest version of the 
Mmtz-Arakawa model (Arakawa, et al., 1974) is based on 
earlier versions developed at UCLA, but contains greater 
vertical resolution extending to the stratopause, an 
improved finite difference formulation, and new parameter- 
lzations of the sub-grid scale processes. While not 
developed for use as a prediction model, numerous existing 
prediction models are based on earlier versions of this 
model and prediction models are being developed which are 
based on. this latest version. What follows is a description 
of the nominal (twelve level) version as of mid-1976. 

MODEL DESCRIPTION 

1.0 Variables 

1.1 Independent 

The independent variables are A and f of a spherical 
polar coordinate system m the horizontal? sigma in the 
vertical defined m the troposphere as 
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where P^. = 100 mb and ? s is the pressure at terrain level, 
and above 100 mb as 


a = 



[M.2] 


where P T = 1 mb; and time, t. 

1 . 2 Dependent 

1.2.1 Prognostic 

The prognostic variables are the pressure, it L , the 
horizontal velocity, V, the temperature, T, and the water 
vapor, q. Secondary prognostic variables are the ground 
temperature, the ground water storage, the mass of snow 
on the ground, the amount of atmospheric ozone, the 
planetary boundary layer depth, and the magnitudes of the 
temperature and moisture discontinuities at the top of the 
boundary layer. 

1.2.2 Diagnostic 

The diagnostic variables are the geopotential , <$ , 
and the vertical velocity, o = . 

2 . 0 Domain 

2.1 Horizontal 

The horizontal domain is a global latitude, longitude 
grid with X and <i> as the coordinates and the grid resolution 
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is 4° m latitude and 5° m longitude. The poles are 
considered as mass points and are treated separately. 

The variables are distributed m the horizontal according 
to Arakawa's Scheme C as shown m Figure M-l. 
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Figure M-l 

Horizontal Grid Structure of the UCLA 12-Level Global Model 


2.2 Vertical 

The vertical coordinate is sigma as defined by Equations 
[M.l] and EM. 2]. Figure M-2 shows the hybrid vertical coordinate 
system of the model consisting of five sigma layers m the 
troposphere and seven constant pressure surfaces m the 
stratosphere. The lowest four layers have equal depth 
m pressure, and the uppermost seven layers have equal depth 
m the logarithm of pressure. 
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Figure M-2 

Vertical Structure of the UCLA 12-Level Global Model 
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2 . 3 Temporal 


Two time levels of data are required since a combination 
of centered and Matsuno time steps is used. A computational 
cycle consists of five six minute time steps with four 
centered steps followed by a Matsuno time step. There is 
no staggering of variables xn time. 

3 . 0 Governing Equations 

The governing partial differential equations written 
m flux form with a sigma coordinate m the vertical are 
listed below. 

Equation of motion: 

ir|| + fkxirV + V(frf) - = irF [M.3] 


Pressure tendency equation: 




V • (rV)da 


First law of thermodynamics : 


[M.4] 


g|; (r^ne) + v* (7rV£n8) + (ira£n0) = 


[M.5] 
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Water vapor equation: 


-- (irq) + V • {irVq) + ^ (nog) = tt(-C+E) 

Equation of state: 
pa = RT 

Hydrostatic, equation: 

ii = ~ RT 

So o 

In the above equations , 

V = horizontal velocity 

f = Coriolis parameter 

k = vertical unit vector 

V = two dimensional gradient operator 

a = specific volume 

"l = P s - P I 

F = horizontal frictional force 
R = gas constant for dry air 

T = temperature 

8 = potential temperature 

= specific heat of dry air at constant pressure 
Q = heating rate per unit mass 


[M. 6] 


[M.7] 


[M.8] 
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$ = gz = geopotential 

q = water vapor mixing ratio 

C = rate of condensation 

E = rate of evaporation 

4 . 0 Numerics 

Finite difference methods are used to approximate the 
solution of the forecast equations. 

4.1 Spacial Differencing 

The distribution of variables on the horizontal grid 
corresponds to Arakawa's (1974) Scheme C as shown m Figure 
M— 1 . The horizontal differencing is constructed so as to 
maintain quadratic constraints analogous to the conservation 
of energy and mean square vorticity. Second-order, centered 
differencing is used m both the horizontal and vertical 
difference operators. 

The change m the variables defined at the poles 
(ir L , T and q) is calculated by summing the meridional flux 
at all points on the latitude circle around the poles. 

4 . 2 Temporal Differencing 

Time differencing proceeds in half-hour cycles using a 
six minute time step. The first four steps of a cycle are 
centered time steps and the fifth is a Matsuno time step. 
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For computational economy, physical processes other than 
advection are only calculated every half hour of simulated 
time on the last time step of a cycle. 

To avoid reducing the time step due to the convergence 
of the grid near the poles, the zonal mass flux and zonal 
pressure force are smoothed longitudinally. A three point 
averaging technique is used for latitudes such that 
3 8 ° 5. I $ I S. 70% and a Fourier filtering technique is used 
for latitudes where j <f) j >_ 74°. No smoothing is performed 
on latitude circles equatorward of 38°. 

A characteristic of the model which is believed to be 
due to numerical truncation occurring in the finite difference 
computations, is a decrease in the globally averaged value 
of surface pressure with time. Therefore, a procedure 
is invoked to maintain the globally averaged value of 
surface pressure at each time step. 

4 . 3 Boundary Conditions 

The horizontal boundary conditions are periodicity m 
longitude and the poles are treated as mass points as 
explained m Section 4.1. 

The vertical boundary conditions are a = 0 at a = -1 
and 1, and the upper pressure P is taken to be 1 mb. 
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5.0 Processes and Effects 


5 . 1 Dynamics 

Terms representing the horizontal and vertical 
advection and convergence/divergence of heat, moisture 
and momentum are included as are terms representing 
the pressure gradient and Coriolis forces m the momentum 
equations. Adiabatic temperature changes are also modeled. 

Orography that is consistent with reality and the grid 
resolution is included to add realism to the flows. 

Surface effects are evaluated as part of the planetary 
boundary layer parameterization which is based on the 
work of Deardorff (Arakawa et al., 1974). 

5.2 Heat Sources/Sinks 

The planetary boundary layer parameterization predicts 
both the height of the boundary layer and the temperature 
discontinuity at the top of the PBL. Quantities m the PBL 
are computed diagnostically. A surface discriminator is 
used to determine which type of surface is being considered 
(land-sea-ice) . The ground temperature is predicted from 
the net surface energy flux. The ocean temperature varies 
spacially, but is held fixed m time. 

Heating due to solar and terrestrial radiation is 
considered including ozone heating m the stratosphere. 
Cloudiness and surface effects are included. 
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Release of latent heat due to large scale precipitation 
is included. 

A detailed parameterization of cumulus ensembles 
(moist convection) which is coupled to the PBL computations 
redistributes heat and moisture m the vertical. A dry 
convective adjustment procedure is included to preclude 
hydrostatic instability. 

5.3 Moisture Sources/Sinks 

The large scale precipitation algorithm removes 
supersaturated states from the model atmosphere and results 
m evaporation in non- saturated layers below the 
saturated layer and/or precipitation. 

The cumulus parameterization process can result m 
the evaporation/condensation of moisture in the lower model 
layers . 

The PBL parameterization predicts both the height of the 
boundary layer and the moisture discontinuity at the top of 
the PBL. Quantities m the PBL are computed diagnostically. 
Evaporation/condensation at the surface is evaluated and 
the soil moisture and snow cover is predicted. 
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6 . 0 Inputs 


Since the UCLA general circulation model is used for 
numerical simulation experiments, and not for weather 
prediction, data is self generated by going through a 
spin-up period from an atmosphere at rest to generate 
"realistic” meteorological conditions. Experiments are then 
made using this initial data. 

7 . 0 Outputs 

An output cycle may be specified at any desired interval. 
Many outputs are built into the code in either grid-point 
format of a field, or printout of various diagnostic 
quantities. In addition to these types of output, the 
entire data base is saved on an output cycle for use 
either m restarting the computation, or for further 
processing. 

8 . 0 Computer Information 

The UCLA general circulation model is programmed 
mostly in Fortran IV (Level H Compiler) for the IBM 360/91 
Computer System. Some specialized output routines are 
coded in assembly language. The execution time was 
approximately 2 hours and 40 minutes per model day for 
the version running in early 1976. 
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N. ANMRC 5-Layer Hemispheric, Spectral Model 

The ANMRC model was developed to study global-scale 
dynamics and numerical weather prediction and is an ex- 
tension of Bourke's (1972) earlier one-level, primitive 
equation spectral model. The version of the model described 
here is dry and adiabatic although there is no fundamental 
limitation precluding the inclusion of detailed physical 
parameterizations. The model also uses a semi-implicit time 
integration of the spectral form of the governing equations 
to obtain a more efficient computation. 

MODEL DESCRIPTION 

1.0 Variables 

1.1 Independent 

The independent variables are the longitude, A. , the 
latitude, <p , sigma defined as: 

» 2- .[N.l] 

P 

* * 

where p is the local surface pressure and time, t. 

* 

1.2 Dependent 
1.2.1 Prognostic 

The prognostic variables are the stream function, i jj , 

the velocity potential, X , the temperature, T and the 

logarithm of surface pressure, q = JLn (p ). 

* 
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1.2.2 Diagnostic 


The diagnostic variables are the scalar velocities 
defined by 

U = ucoscf) and V = vcos<f> [N.2] 

the geopotential , $ , and the vertical velocity, a. 

2 . 0 Domain 

2.1 Horizontal 

The horizontal domain is hemispheric (or global) with a 
J=15 wavenumber truncation (rhomboidal) coupled with a 
hemispheric transform grid of 20 Gaussian latitudes and 48 
equispaced longitudes. The poles are regular points m this 
representation and variables are placed regularly on the 
transform grid. 

2 . 2 Vertical 

The vertical coordinate is sigma as defined by Equation 
[N.l] . Figure N-l shows the model sigma surfaces and the 
placement of variables in the vertical. 
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Figure N-l 

Vertical Structure of the ANMRC 5-Layer Spectral 
Model (Repeated for all levels) . 


2 . 3 Temporal 

Two time levels of data are required (except during the 
startup step) since a centered semi -imp licit time differencing 
scheme is used with a sixty minute time step. There is no 
staggering of the variables m time. 
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3.0 Governing Equations 


The governing partial differential equations written m 
spherical polar coordinates m the horizontal and a sigma 
coordinate m the vertical are listed below: 

Vorticity equation: 


9 

9 t 


V 2 i|> 


1 ,9 A , , 9B, 

( 9T + COS(t> 


2 fl (sin|V 2 X + J) 

d 


+ K V 2 (V 2 ip 

li 


2$ 

— 2" 


[N.3] 


Divergence equation: 


9 1 


“ v 2 x 


9B 


acos^ '3X 


- 003*11) + 


2 q (sm^V 2 ^ - 


1) 


9 T 
9 1 


9q 

9 1 


- V 2 (E + $’ + RT q) 

o 


Thermodynamic energy equation: 


(2S£- + ”*1) + T ’ V 2 x 


, acos cp 9 A 


~w 


RT 


+Y0 + ~ tV 2 X + (V + V)*Vq] + K h V 2 T 


Pressure tendency equation: 


V-Vq + V~ 2 x + K 0 V 2 q 
H 


[N.4] 


[N.5] 


[N.6] 


N-4 



Hydrostatic relation: 


3$ RT 

3a a 


[N.7] 


Here : 


U 


cos(p 3 i|j 1 3X 

a 3 a 3 X 


[N.8] 


= 1 li + cos<j> 9X 
a 3X a 3<J) 


[N.9] 


kxVif> + VX 


(l-a)jf V 2 Xda ~ J^Vxda 


+ [ (1-a) £ Wda - j^Vda]*Vq 

[N. 10} 

, * 3V , RT" , 3q q 3X y 

V * V + °3o + a 3* - p* 3o 

[N.11I 

v**v - i su RT ' 8q + ~^ 9Tx 

v vv a 3a a 3X P* 3a 

[N.12] 

2 2 
u z + v z 

2 . 

[N.13] 


_ 2 
2 cos <j) 


the subscripted zero denotes a horizontal mean value and the 
superscripted prime the deviation from that mean; and the 

r 0 

bar notation means ( ) = { ) da . 

In the above equations, 

U = ucos<|> = scalar zonal wind component 

V = vcos<J> = scalar meridional wind component 

= stream function 
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velocity potential 

£n{p*) where p* is the surface pressure 
T q + T ' , temperature 
vertical velocity 

= static stability 

crCp Sa 

+ $* = geopotential height 
rotation rate of the earth 
horizontal stress m the x direction 
horizontal stress m the y direction 

4.0 Numerics 

A spectral representation m the horizontal coupled with 
finite differencing m the vertical and semi-implicit time 
differencing is used to approximate the solution of the model 
equations. 

4.1 Spacial Representation . 

The essence of the transform procedures proposed by 
Orszag and Eliasen is to evaluate the nonlinear terms of 
the equations as simple products following a transform from 
the spectral to gridpomt domain; a subsequent inverse 
transform yields the spectral form of the requisite nonlinear 
term. The procedure may be more generally described as com- 
bining nonlocal (spectral) differentiations with local 
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(gridpomt) multiplications for evaluation of the nonlinear 
terms . 

That is, the model considers the prognostic variables 
% X and T m terms of orthogonal spherical harmonics at 
each of five levels m the vertical; the diagnostic varia- 
bles U and V are similarly represented at these levels 
as is the prognostic variable q at the surface. Truncated 
expansions for the prognostic and diagnostic variables are 
required as follows ; 


ip - 


J |m |+J 
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m=-J £=|m[ 
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JL l 


[N.14] 


where 

a) The terms iJj™, etc. denote time dependent, generally 
complex, expansion coefficients, and the reality of the fields 
requires = (-l) m and so forth. 
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b) The spherical harmonics are given by Y™ = P™ (sm<|>} e im ^ 
where P™ (sintft) is an associated Legendre polynomial of 

the first kind normalized to unity. 

c) The wavenumber J denotes rhomboidal truncation. 

Diagnostic relationships for U™ and V 1 ^ as obtained from 
Equations [N.8] and [N.9] are: 

= (H)eJ -(£+2)e™ +1 ^ +I + im X® [N.15] 

V 2 = - + U*2)e m l+1 X“ +1 + WJ tH.16] 

where 

e™ = IU 2 - m 2 ) / (4 (£) 2 - 1} ] 1/2 

The spectral representation of nonlinear products 
occurring m the prognostic and diagnostic equations requires 
spectral grid transforms. The nonlinear products are evaluated 
on a latitude-longitude grid (equally spaced m longitude, 
Gaussian in latitude) at each model level. These products 
are then expressed as truncated Fourier series at each 
latitude circle giving: 
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V*4,u + ^ * 5|1 cos*|a - 2- ^ = a jr 
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[N. 17 ] 


VT* = a 
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'V 2 X + <?Y + ^ V 2 xda + (V + Vdcr)*Vq] = 
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m=- J 
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xmA 


m=-J 


The dxagnostxc evaluatxon of a xs 'requxred only at gridpoxnts 
and xs not subsequently expressed spectrally. 
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With additional representations given by 

( = J/ ( ) m pJ (sin*) cos* d * [N. 18] 

~2 

where the integral is exactly evaluated by Gaussian quadra- 
ture and applied to E m , and Z m thereby yielding E^, 

and Z™, the equations for the spectral amplitude ten- 
dencies for the stream function, velocity potential, tem- 
perature and logarithm of the surface pressure are found 
to be : 


- 1 ( 1 + 2 .) 
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[N. 21] 
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m 

■ z 7- £<i+1) / 1 x " ac ’ 

The vertical derivatives appearing m the various non- 
linear terms and the vertical integrals are all evaluated by 
finite differencing at the half-levels. Full-level values are 
obtained by averaging half-level values where appropriate. 
Logarithmic variation of temperature between the lowest two 
levels is assumed m defining at the lowest full and 
half- levels. 

4.. 2 Temporal Differencing 

A semi-implicit time integration algorithm based on the 
worlc of Robert (1969") is used m the model. The basic idea 
of the semi-implicit technique is to treat the gravity wave 
modes implicity (unconditional stability) to give a much 
longer time step. This results in a time step of 60 minutes. 

In order to facilitate description of the semi-implicit 
algorithm, it is convenient to rewrite the spectral prog- 
nostic equations in a simpler form. The prognostic equations 
for the velocity potential, temperature and logarithm of the 
surface pressure may be written as r 

3D™ 

+ IU+1) + R'T^q™) [N. 23] 
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[N.2M 
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= y z + 


[fof da 
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+ JT° o”da 


[N. 25] 


where D™ = -£(£+1) X^, the quantities x™, y™ and z™ denote 

the spectral representation of the nonlinear terms and other 

2 

terms not explicit here, and R' denotes R/a . 

By writing the equations m matrix form where the 
vectors are column vectors over the number of levels, the 
prognostic - equations may be written as : 


si 

3 1 


= R'A y + AcrR' A G D 


[N. 26 ] 


3D ' 

^rr = X + ZU+ 1) [$ + R'T q] - [N. 27] 

o'- — — — o 


If = z + ActI D \ [N. 28] 

where $_ = R'A T has been used to obtain the first equation. 

Using a centered time difference and the semi-implicit 
operator of Robert, the finite difference equations may now be 
written as : 

T+l T— 1 T+l T— 1 

£ = i + 2AtR' A y; + AtAcrR' AG (D + D ) [N.29] 
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D x+1 = D t 1 + 2 Atx + At 1(1+ 1) [(£ T+1 + $ T 1 ) 

+ R'T q (q T+1 + q^ 1 )] [N.30] 

q T+1 = q T_1 + 2Atz + AtAal (D T+1 + D^ 1 ) [N.31] 

These three equations involve three unknowns that may be 
solved for the variable D. Then by using Gauss elimination, 

' the system is solved for D and the other variables obtained 
by back substitution. 

A forward time step is performed on the first time step. 

4 . 3 Boundary Conditions 

The representation of the variables in terms of spherical 
harmonics implies periodicity m longitude, and the poles, 
although not transform grid points, are regular points. 

The vertical boundary conditions are a = 0 at o = 1 
and 0. 

5.0 Processes and Effects 

5 . 1 Dynamics 

Terms representing the horizontal and vertical advection 
and convergence/divergence of heat, vorticity and divergence 
are included as are terms representing the pressure gradient 
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and Coriolis forces m the vorticity and divergence equations. 
Adiabatic temperature changes are also modeled. 

Orography that is consistent with reality and the 
transform grid resolution is included to add realism to the 
f lows . 

Horizontal stress is applied at the lowest level to 
simulate vertical momentum transport at the surface. 

5.2 Heat Sources/Sinks 

A dry, adiabatic fluid is considered? thus, there are 
only adiabatic heating terms. 

5.3 Moisture Sources/Sinks 

There are none since a dry, adiabatic fluid is considered. 
6 . 0 Inputs 

The model was originially tested with Phillips analytic 
initialization and found to have good phase speed charac- 
teristics as was expected of the spectral representation. 

The global integrals of total energy, angular momentum and 
mass showed good conservation properties. 

Test integrations were performed with GARP data sets 
from the Melbourne WMC as analyzed by the method of Gauntlett, 
et al. (1972) which gives temperatures and winds at standard 
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levels and the surface pressure. The analyzed data is then 
interpolated to the transform lat-long grid and interpolated 
to the sigma surfaces. Fourier and associated Legendre 
transforms are then performed to obtain the initial spectral 
representation of the data fields. 

7 . 0 Outputs 

Fields at standard levels and the surface pressure are 
produced at variable output times. 

8 . 0 Computer Information 

The model is coded m FORTRAN (level-H Compiler) for 
the IBM 360/65 (185K bytes of storage required) . The 
5- layer, J=15 model requires 13 minutes to perform a 24-hour 
hemispheric forecast. 
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O. Canadian 5-Layer Hemispheric Spectral Operational Model 

The model described here became the current operational 
large-scale forecast model of the Atmospheric Environment 
Service of Canada on 18 February 1976. The model traces 
its development to the work of Orszag and Eliassen (c. 1970), 
Bourke (1974) of the ANMRC and Robert (1972) of Canada 
and is probably the most complete spectral primitive 
equation model in existence today. 

MODEL DESCRIPTION 

1.0 Variables 

1 . 1 Independent 

The independent variables are the longitude, X, the 
latitude, <j>, sigma defined as 



where p* is the local surface pressure and time, t. 

1 . 2 Dependent 
1.2.1 Prognostic 

The prognostic variables are the vorticity , £, the 
divergence, D, an energy quantity, P = $ + RT^q, the 
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logarithm of surface pressure, q = £n(p*) , and the dew 
point depression, S = T - T^. A secondary prognostic 
variable is the precipitation. 

1.2.2 Diagnostic 

The diagnostic variables are the scalar velocities 



V 


V cos (j) 
a 


[ 0 . 2 ] 


the geopotential , $ , and the vertical velocity, a. 

2 . 0 Domain 

2.1 Horizontal 

The horizontal domain is hemispheric with a J = 20 

wavenumber truncation (rhomboidal) coupled with a transform 

grid of 26 Gaussian latitudes and 64 equispaced longitudes. 

* 

The poles are regular points m this representation and 
variables are placed regularly on the transform grid. 

2 . 2 Vertical 

The vertical coordinate is sigma as defined by 
Equation [0-1J * Figure 0-1 shows the model sigma surfaces 
and the placement of variables m the vertical. 
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Figure 0-1 

Vertical Structure of the Canadian 
5-Layer Spectral Model (Repeated for all levels) 
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2.3 Temporal 


Two time levels of data are required (except during the 
startup step) since a centered semi-implicit time differencing 
scheme is used with a forty minute time step. There is no 
staggering of the variables m time. 

3 . 0 Governing Equations 

The governing partial differential equations written 
m spherical polar coordinates in the horizontal and a 
sigma coordinate in the vertical are listed below. 

Vorticity equation: 


|| = -a(A,B> 
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Divergence equation: 


3D 
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Energy equation: 
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Pressure tendency equation: 


H - w s - 0 


[ 0 . 6 ] 
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Moisture equation: 


|| = -a (US , VS) + B s 


Hydrostatic relation: 
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[0.17] 
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and the subscripted zero denotes a horizontal mean value 
and the superscripted prime the deviation from that mean. 
In the above equations , 

f = Coriolis parameter 

„ u cos<f> , , , 

U = — — — = scalar zonal wind component 

Y _ v -^ - ° - s ^ = scalar meridional wind component 

z> = vertical component of vorticity = k*VxV 

D = horizontal divergence = V*V 

T = T + = absolute temperature 

o 

q = £n(p*) = logarithmic pressure 

Y = Y q + Y'" = ~~ - = static stability 


o 




= vertical motion m sigma coordinates 
= gz = geopotential 

= horizontal components of frictional force per unit mass 
= diabatic heating 
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R = gas constant for dry air 

0^ = the specific heat of dry air at constant pressure 

$ 

= dewpoint temperature 

S = t - is the dewpoint depression 

e = ratio of the molecular weight of water vapor to 

effective molecular weight of dry air (0.622) 

LfT^) = latent heat of vaporization of water and ice 

H.. = moisture sources or sinks 

M 

4.0 Numerics 

A spectral representation m the horizontal coupled 
with finite differencing m the vertical and semi-implicit 
time differencing is used to approximate the solution of 
the model equations. 

4 . 1 Spacial Representation 

The essence of the transform procedures proposed by 
Orszag and Eliasen is to evaluate the nonlinear terms of 
the equations as simple products following a transform from 
the spectral to gridpomt domain; a subsequent inverse 
transform yields the spectral form of the requisite non- 
linear term. The procedure may be more generally described 
as combining nonlocal (spectral) differentiations with local 
(gridpomt) multiplications for evaluation of the nonlinear 
terms . 
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The spectral form of the governing equations is developed 
m a manner analogous to that of Machenhauer and Daley (1972) 
and, more specifically, Bourke (1974). The vorticity, §, 
for example, is expanded in the following truncated series 
of spherical harmonics: 

J | m 1 + J 

C = 2 2 ^ [°-19] 

m=-J l=\m\ c 

where 


a) the terms denote time dependent, generally complex, 
expansion coefficients, and the reality of the fields 
requires ( C J) * = (-l) m Z~ m 


b) the spherical harmonics are given by = P™ (sm<j>)e 
where P™ (sm<f>) is an associated Legendre polynomial of 
the first kind normalized to unity 

c) the wavenumber J denotes rhomboidal truncation. D, T, 
P , $ , S and W are expanded in the same manner as Equation 
[0.19] . The variable q is similarly expanded except that 
the q^ are functions of time only. 

In the expansions for U and V there is one extra 
component for each m, m order to be consistent with the 
expansions for C and D (Eliasen et al. 1970) . 

J | m | + J+l J |m|+J+l 


o-2 1 

m=-J £=jm| 


n m v m 
l *1 


V = 



m=-J 



V* Y™ [0.20] 
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Diagnostic relations between U, V and D, £ can be obtained 
through Helmholtz's theorem which states that an arbitrary 
wind field can be represented as the sum of the gradient 
of a potential function and the curl of a stream function; 
thus 

IU+ l)uj? = -U+ 1 )eJe 2 _ 1 + -te“+l c £+l " imD l [0.21] 

and 

IU + DVJ = (£+1)eJdJ_ 1 - £e^ +1 D™ +1 “ EO- 22] 


where 


el = E (£ 2 ~ m 2 )/( U 2 - 1)] 1/2 Eo.23] 

U° = e°C° and V° = -e°D° are special cases, 
oil o 11 

Other diagnostic relations can be obtained from 
the hydrostatic equation; thus 
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and 


T m 
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[0.25] 


0-9 



The prognostic variables are z™, , P™, S™ and q™ . 

The variables T^, and V™ are all related to the 

prognostic variables by the diagnostic relations. 
Introducing the operator 






(A,(j>) cos<i>dAd<}> 


[0.26] 


the spectral form of the governing equations become 
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[0.27] 

[0.28] 

[0.29] 

[0.30] 

[0.31] 

[0.32] 


All variables required in the calculation of the 
right-hand sides of Equations [0.27-0.32] are first 
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synthesized onto the physical space transform grid. The 
nonlinear expressions are then evaluated and the integrals 
on the right-hand side are then simply calculated by exact 
numerical quadrature. 

The transform grid consists of equally spaced longitudes 
and almost equally spaced Gaussian latitudes. In the 
synthesis from spectral to grid, an associated Legendre 
transform is first performed at each Gaussian latitude 
circle to produce Fourier coefficients. A discrete 
Fourier transform is then performed to produce functional 
values at each longitude of the transform grid. After the 
nonlinear calculations are carried out on the grid, the 
integrals on the right-hand sides of the equations are 
calculated by first performing a discrete inverse Fourier 
transform at each Gaussian latitude. The final step is 
the inverse associated Legendre transform performed by 
Gaussian quadrature. 

The vertical finite-difference scheme is a somewhat 
more general form of the scheme used by Robert et al . (1972). 
The basic feature of the scheme is that the temperatures 
are carried at levels intermediate to the levels of the 
geopotentials . The vertical finite differencing involved 
in the right-hand sides of Equations [0.27-0.32] is 
conventional except that logarithmic differencing has 
generally been used, particularly with respect to the 
thermodynamic variables. 
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4.2 Temporal Differencing 

A semi-implicit time integration algorithm following 
Robert et al. (1972) is used m the forecast model. The 
basic idea of the semi-implicit technique is to treat the 
gravity wave modes implicitly (unconditional stability) 
and the meteorological modes explicitly (conditional 
stability) to give a much longer time step. This results 
m a time step of forty minutes. 

The time derivative is approximated by a centered 
difference (except at the initial time where a forward 
difference is used) . The remaining terms on the left-hand 
sides of Equations [0.27-0.32] are handled implicitly by 
the application of the time averaging operator defined by 

¥ t = 1/2 [F (t + At) + F(t - At)] [0.33] 

while the right-hand sides are calculated explicitly. 

Two equations remain fully explicit. They are the 
vorticity and moisture equations: 

eJ(t+At) = -2At{a(A,B) + ej (t-At) [0.34] 

and 

sJ(t+At) = 2At{-a (US , VS) + B + S™ (t-At) [0.35] 
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The remaining equations are handled m a semi-implicit 


manner giving: 



At t (£+1) 
a 2 


pJ = At{a(B,-A) 
+ D™ (t-At) 


a 2 V 2 E}^ 


gp'" t 

a L- + AtRY o Wj = AtR{ct(UT%VT") - 

3 p » 

+ a IhT 

\ 

— t 1 

m 4I _ „ in m A»\ 

q l ■ At W s l = q l 

and 


[ 0 . 36 ] 


[ 0 . 37 ] 


[ 0 . 38 ] 


3W 


m 


3a 


+ D 


m 


<V" 


[ 0 . 39 ] 


The integrals on the right-hand side are all 
evaluated at time t. Some of the terms m the integrals 
are forcing terms and should really be calculated at time 
t-At to avoid the growth of a computational mode. This 
would require considerable additional computation, so 
the forcing terms are all calculated at time t and the 
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computational mode is suppressed with a very weak time 
filter (Asselm, 1972) . 

The implicit system of equations [0.36-0.39] is 

— t — — t 

solved by eliminating D™ ? then is eliminated from 

the resulting equations giving a diagnostic relation for 
— t 

P™ . This equation can then be solved for P™ and the 
remaining variables obtained by back substitution. 

4 . 3 Boundary Conditions 

The representation of the variables m terms of 
spherical harmonics implies periodicity in longitude, 
and the poles, although not transform grid points, are 
regular points . 

The vertical boundary conditions are a = 0 at a = 1 

and 0 . 

5.0 Processes and Effects 

5.1 Dynamics 

Terms representing the horizontal and vertical 
advection and convergence/divergence of heat, moisture, 
vorticity and divergence are included as are terms 
representing the pressure gradient and Coriolis forces m 
the vorticity and divergence equations. Adiabatic 
temperature changes are also modeled. 
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Orography that is consistent with reality and the 
transform grid resolution is included to add realism to 
the flows . 

The skin friction F and F, are calculated at the 

A <p 

lowest level only using Cressman's drag formulation. 

5.2 Heat Sources/Smks 

Sensible heating m the PBL is evaluated at the 
lowest model level over the open ocean. The term 
depends on the wmdspeed and the difference between 
the extrapolated surface temperature obtained from the 
lowest sigma levels and the sea surface temperature. 

The sea surface temperature is spacially varying, but 
held constant during the forecast. 

Release of latent heat due - to large-scale precipi- 
tation is included. 

A moist convective adjustment process based on the 
method of Manabe et al. modifies unstable vertical profiles 
of temperature and moisture conserving moist static 
energy. A dry convective adjustment procedure is included 
to preclude hydrostatic instability. 

5.3 Moisture Sources/Smks 

The large scale precipitation process removes super- 
saturated states from the model atmosphere. 
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A moist convective adjustment process based on the 
method of Manabe et al . modifies unstable vertical profiles 
of temperature and moisture conserving moist static energy. 
This process may result m convective precipitation. 

Evaporation in the PBL is evaluated at the lowest 
model level over the open ocean. The term H M depends 
on the wmdspeed and the difference between the lowest 
layer mixing ratio and the saturated mixing ratio of air 
m contact with the sea. 

6 . 0 Inputs 

The forecast model is initialized from an objective 
analysis (Rutherford, 1976) which provides winds, heights, 
temperatures and dewpoint depressions at standard pressure 
levels and surface pressure on a 2805 point polar stereo- 
graphic grid (381 km grid length at 60°N) . This data is 
then interpolated with bi-cubic splines to the Gaussian 
transform grid. Some extrapolation is necessary m the 
tropics, and the symmetry conditions are used to determine 
the Southern Hemisphere. 

The surface topography is then used to vertically 
interpolate the data to the model sigma suf faces. Finally, 
the spectral coefficients of the model variables are 
obtained by spectral analysis of the sigma surface data. 
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7 . 0 Outputs 


The sea level pressure, precipitation and heights 
and temperatures at standard pressure levels are produced 
on a normal output cycle . 

8 . 0 Computer Information 

The forecast model is programmed mostly in CDC 
Fortran Extended (a dialect of FTN IV) for the CDC 7600. 
Approximately 25k words of SCM and lQOk words of LCM 
are required for the SP20 hemispheric model. The 
running time is approximately five minutes for a 36 hour 
forecast. 


9.0 References 


Asselm, R. , "Frequency Filter for Time Integrations," 

Monthly Weather Review, Vol. 100, No. 6, June 1972, 
pp. 487-490. 

Bourke, W. , "A Multi-Level Spectral Model. I - Formulation 
and Hemispheric Integrations , " Monthly Weather 
Review, Vol. 102, No. 10, October 1974, pp. 687-701. 

Daley, R. , C. Girard, J. Henderson and I. Simmonds , 

"Short-Term Forecasting with a Multi-Level Spectral 
Primitive Equation Model," Atmosphere, Vol. 14, 

No. 2, 1976, pp. 98-134. 

Eliasen, E.B., B. Machenhauer and E. Rasmussen, "On a 

Numerical Method for Integration of the Hydrodynamical 
Equations with a Spectral Representation of the 
Horizontal Fields," Institute of Theoretical 
Meteorology, University of Copenhagen, Report No. 2, 
1970. 


0-17 



Robert, A. , J. Henderson and C. Turnbull, "An Implicit 
Time Integration Scheme for Baroclinic Models 
of the Atmosphere," Monthly Weather Review, 

Vol. 100, No. 5, May 1972, pp. 329-335. 

Rutherford, I.D. , "An Operational Three-Dimensional 
Multivariate Statistical Objective Analysis 
Scheme," JOC Study Group Conference on Four- 
Dimensional Data Assimilation - Paris, 17-21 
November 1975, Report No. 11, the GARP Programme 
on Numerical Experimentation, 1976, pp. 98-121. 


0-18 



P. Deutscher Wetterdienst 6-Layer Hemispheric Operational Model 


The Deutscher Wetterdienst atmospheric prediction 
model is currently executed twice daily to produce four 
day forecasts for the Northern Hemisphere. Although 
the operational model described here is a dry model, 
a research version of the model also exists which 
contains moisture. The following description represents 
a snapshot of the state of the operational model as of 
mid-1976. 

MODEL DESCRIPTION 

1 . 0 Variables 

1 . 1 Independent 

The independent variables are x and y derived from a 
polar stereographic map projection in the horizontal, 
pressure m the vertical, p, and time, t. 

1 . 2 Dependent 
1.2.1 Prognostic 

The prognostic variables are the x wind component, u; 
the y wind component, v; and, the potential temperature, 0. 
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Secondary prognostic variables are the geopotential of 

the 1000 mb surface, ^iq 00 ? anc ^ t * ie surface temperature 

over land or ice, T . 

s 

1.2.2 Diagnostic 

The diagnostic variables are the geopotential , $ , 
and the vertical velocity m pressure coordinates, w. 

2 . 0 Domain 

2.1 Horizontal 

The horizontal domain is a polar stereographic projection 
(true at 60°N) of the Northern Hemisphere north of 12 °N 
(standard octagon) using x - y coordinates with map factor 


1 + sin 60° 
m ~ 1 + sm $ 


[P.l] 


The grid resolution is 381 km at 60°N and the grid is 
shown in Figure P-1. 

The variables are staggered m space and time according 
to the scheme of Eliassen (1956) which is shown m Figure P-2. 
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Figure P-1 

Northern Hemisphere Polar Stereographic Grid 
(standard Octagon) 


P-3 




even timestep 


odd timestep 


Figure P-2 

Horizontal Staggering of Variables in Space and Tune 
for the Deutscher Wetterdienst Model 
(Eliassen Grid) 





2 . 2 Vertical 


The vertical coordinate is pressure. Figure P-3 
shows the vertical structure of the model and the 
placement of variables m the vertical. 
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Figure P-3 

Vertical Structure of the Deutscher Wetterdienst Prediction 

Model 


2 . 3 Temporal 

Two time levels of data are required since a centered 
time differencing scheme is used with a five minute time 
step. The variables are staggered m space and time according 
to the scheme of Eliassen (1956) as shown in Figure P-2. 
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3 . 0 Governing Equations 


The governing partial differential equations, written 
m Cartesian coordinates in the horizontal and a pressure 
coordinate m the vertical, are listed below. 

Wind component m the x direction: 

!“ = - U |H _ V |S _ + fv - M + D + F IP.; 

9t 3x 3y 3 p 3x u x 


Wind component m the y direction: 


- U |V - v|X - JZ - f u - |» + D + F 
3x 3y 3P 3y v y 


Thermodynamic energy equation: 


,,96 ..9 0 3 0 _ , rj 

' " v 8y ‘ “3? + D T + H 


Continuity equation: 


V*V + |^- = ° 
3p 


Hydrostatic equation: 


Equation of state: 


p = pRT 


[P.7] 
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D m = 


above equations, 
x-component of wind velocity 
y-component of wind velocity 
horizontal velocity vector 
vertical (p) velocity 
potential temperature 
gz — geopotential 
Coriolis parameter 
horizontal diffusion terms 

stress components 

sensible heating term 
density 

gas constant for air 


4.0 Numerics 


Finite difference methods are used to approximate the 
solution of the forecast model primitive equations . 

4.1 Special Differencing 

Centered, second-order finite difference methods are 
used to approximate the horizontal and vertical derivatives. 
Since the staggered grid system of Eliassen is used, 
averaged values of the variables are required m the finite 
difference approximations. As an example of this procedure/ 
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the finite difference form of the thermodynamic energy 
equation is given by 


90 — *y _ x _ v _ v ■■■ - P — xy 

It = " m 2 [ U X 0 Y + V Y 0 X + Q0 Xp ] + m 2 A h V 2 0 [P.8] 

where U = — , V = ^-r , ft = — and the bar notation 
m 2 m 2 m 2 

denotes averaging in the direction specified. 

, The lateral diffusion terms in the x and y component 
wind equations and thermodynamic energy equations take 
the form of 

Di|> = AV 2 i|> [P.9] 

where both the n-1 and n time levels are used in the 
computation of the laplacian. 

4 . 2 Temporal Differencing 

Centered, second-order time differencing is used to 
approximate the time derivatives with a time step of five 
minutes. The computation is started with a forward time 
step. 

4 . 3 Boundary Conditions 

The horizontal boundary conditions are that there is 
no flow normal to the boundary and no heat flux across the 
boundary. 
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The vertical boundary conditions are m = 0 at p = 0 

at the model top, and w = V ‘VZ at the model surface. 
cr s s s 

5.0 Processes and Effects 

5.1 Dynamics 

Terms representing the three-dimensional advectxon of 
momentum and heat are included as are terms representing 
the pressure gradient and Coriolis forces m the equations 
of motion. 

Orography that is consistent with reality and the 
grid resolution is included to add realism to the forecasts. 

Friction based on the vertical derivative of the 
stress is applied at the lowest two model levels. The 
stress at the 1000 mb level is given by a drag coefficient 
formulation and at the 800 mb level by the wind shear. 

5.2 Heat Sources/Sinks 

Sensible heat exchange between the surface (ocean 
as well as continents) and the lowest model level is 
included. The ocean temperature is spacially varying, 
but held constant during a forecast. 

5.3 Moisture Sources/Sinks 

There are none since this is a dry model. 
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6 . 0 Inputs 


The forecast model initial data is obtained from 
hemispheric analyses of the mass structure at the 850, 

700, 500, 300, 200 and 100 mb levels using the iterative 
correction method. In addition, analyses are made of 
the surface pressure (smoothed m the tropics) , the 
temperature at the two meter level and the ocean temperature. 

Winds are obtained from the solution of the balance 
equation on the model pressure surfaces. 

7 . 0 Outputs 

A normal output cycle produces the surface pressure; 
w on the 150, 250, 400, 600 and 800 mb surfaces; and the 
winds, geopotentials and temperatures on the 100, 200, 

300, 500, 700 and 850 mb surfaces. An output cycle is 
made at the initial time plus 18, 24, 30 and 48 hours. 

8 . 0 Computer Inf ormation 

The forecast model is run on a Control Data Corporation 
3800 Computer System. The model is normally executed 
twice daily and requires about two and one-half hours to 
produce a forty-eight hour forecast. 

The model is written mainly in Fortran with some 
routines in assembly language and requires approximately 
32K of central memory and 40 OK of mass storage. 
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Q. Japan Meteorological Agency 6-Layer Limited Area 
Operational Model 

The Japan Meteorological Agency operational fme-mesh 
limited area model (FLM) is executed twice daily to produce 
24-hour prognostic charts for Japan and the surrounding 
area. As with all operational models, changes are contin- 
ually being made to improve the quality of the forecasts 
produced and to better utilize the available computer re- 
sources; thus, this description is a snapshot of the fore- 
cast model as defined by the available reference documenta- 
tion (1974) . 

MODEL DESCRIPTION 

1.0 Variables 

1.1 Independent 

The independent variables are x and y derived from 
a polar stereographic map projection m the horizontal; 
pressure, p, m the vertical; and time, t. 

1.2 Dependent 
1.2.1 Prognostic 

The prognostic variables are the surface pressure, p^, 
the temperature, T, the specific humidity, q, the x wind com- 
ponent, u, and the y wind component, v. A secondary prog- 
nostic variable is the precipitation, P r , obtained from both 
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the large-scale precipitation algorithm and the moist convec- 
tive adjustment (cumulus parameterization) . 

1.2.2 Diagnostic 

The diagnostic variables are the geopotential, <J>, and 
the pressure coordinate vertical velocity, w. 


2.0 Domain 

2.1 Horizontal 

The horizontal domain is a polar stereographic projec- 
tion of an area approximately 5600 by 5000 km centering on 
the Japan Islands using x-y coordinates with map factor: 


1 + sm60 ° 
1 + sm 4> 


[Q.l] 


The horizontal grid resolution is 152.4 km at 60 °N. Figure 
Q-l shows the fine-mesh limited area (FLM) grid nested m the 
coarse mesh grid from which boundary conditions for the FLM 
are obtained. There is no staggering of variables m the 
horizontal . 


2.2 Vertical 


The vertical coordinate is pressure. Figure Q-2 shows 
the model computational levels and the placement of the varia- 
bles m the vertical. 
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Figure Q-l 


Japan Meteorological Agency Fine Mesh 
Area Grid Nested m the Coarse Mesh 
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Figure Q-2 

Vertical Structure of the Japan Meteorological 
Agency Forecast Model 













2 . 3 Temporal 


One time level of data is required since the Euler- 
backward scheme is used. There is no staggering of variables 
in time . 


3 . 0 Governing Equations 

The governing partial differential equations, written 
for an x-y map projection with map factor m, are listed 
below. 

Equation of motion m the x-direction: 


3u 

3t 



[— (— ) 
L 3x ’ 


3_ 

ay 


(^)I 

m 


S (urn) 5 $ 

3p m 3x 


+ fv + m v [v 


8m 
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3m, 
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+ F x + D 


[Q.2] 


Equation of motion m the y-direction: 


3v „2 r 3 ,uv, , 3 ,vv,, 3(vw) „ 3$ 

at “ [ a5 'i -1 + ay <-H> ] 5p~ " m ay 


, j; 2 r 3m 3m, , 

+ fu - m u [v ^ - v -g^] + F y + D 


[Q.3] 


Thermodynamic energy equation: 


3T 2 r 3 ,uT, , 3 ,vT,, 3 (Tic) 

at = “ m [ 3iF ( i“ ) + 37 ( ST n 3i~ 


o 3# 

C 3p 
P F 


Q 


[Q.4] 
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Moisture equation : 


9q _ 

3t 


3P S 

at 


3(i) _ 
3P 


- m 2 [ 3 (HS) + 4_ ( za n - l<Sii + 

3x v m 3y v m ' J 3p 


M 


Surface pressure tendency equation: 


3P S 3P S 

“ m [U + v— 9y^ + w s 


Continuity equation: 


- m 2 (-) + (-) ] 

3x in 7 3y 'm ; 1 


Hydrostatic equation: 


9$ RT 

3P P 

Equation of state : 
p = pRT 
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In the above equations, 

u = x-direction wind component 

v = y-direction wind component 

T = temperature 

q = specific humidity 

to = vertical velocity m p coordinates 

f = Coriolis parameter 
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stress components 


diffusion term 

dxabatic heating/cooling term 

gz = geopotentxal 

densxty of axr 

surface pressure 

gas constant for dry axr 

specxfic heat of dry axr at constant 
pressure 


Fxnxte dxfference methods are used to approxxmate the 
solutxon of the forecast model equations. 

4.1 Spacxal Dxfferencxng 

Centered, quasi-energy conservxng, second-order fxnxte 
dxfference operators are used to approxxmate the horxzontal 
derxvatxves. Vertical derxvatxves are approxxmated by cen- 
tered, second-order dxfferences. 

A horxzontal smoothxng operator xs applxed to the tem- 
perature, geopotentxal, specif xc huraxdxty and surface pres- 
sure frelds to control the computatxonal noise. For the phy- 
sical quantity A, the operator takes the form: 


F , F 
x y 

D 

Q 

$ 

P 

p s 

R 

C 


4.0 Numerics 



5 2 

Where K is the smoothing coefficient (10 xd ) . K is 
increased by a factor of eight and four, respectively, on 
the outermost and next to outermost line from the grid 
boundary. 

The dissipation term m the equation of motion is given 

by: 


D = k-^kxV? + k 2 VV EQ.ll] 

t 

Where C and V are the relative vorticity and horizontal di- 
vergence, respectively. The diffusion coefficients have the 
values k^ = 10 5 and k 2 - 10^ and are increased near the 
boundary in the same manner as the smoothing operator dif- 
fusion coefficient. 

4.2 Temporal Differencing 

Euler-backward time differencing is used to integrate, 
the forecast model equations. This differencing technique 
gives additional control of the computational noise since 
it is dissipative. 

The diabatic heating terms are evaluated every third 
time step to reduce the amount of computation. 

The area mean surface pressure is kept constant through- 
out the forecast integration by adjusting the individual 
surface pressures. 
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4 . 3 Boundary Conditions 


The horizontal boundary conditions for the FLM are 
obtained from the coarse-mesh model (QGM) . A twenty-four 
hour integration is performed with the QGM and the coarse- 
mesh values are saved every hour for use m the FLM inte- 
gration. The FLM boundary values are updated hourly and 
are obtained from the space interpolated QGM values through 
the following procedure: 

a) 4, q, p s and the tangential component of the wind 
with respect to the FLM grid boundary are replaced by the 
QGM values . 

b} T is replaced as m a) at an inflow boundary point 
and is extrapolated from the FLM grid at an outflow boundary 
point. 

c) The normal component of the wind with respect to 
the FLM grid is determined through an application of the 
continuity equation (Equation Q.7). 

In addition, values of p and T m the FLM next to the 

s 

boundary are replaced by the QGM values hourly. 

The vertical boundary conditions are to = 0 at p = 0 

and to = to at the surface, 
s 
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5.0 Processes and Effects 


5 . X Dynamics 

Terms representing the three-dimensional flux diver- 
gence of heat, moisture and momentum are included, as are 
terms representing the pressure gradient and Coriolis forces 
m the equation of motion. Adiabatic temperature changes 
are also modeled. 

Orographic effects are also included to add realism 
to the forecast features. 

Friction evaluated according to the vertical deri- 
vative of the stress is included at the lowest level. The 
stress at the surface depends on a drag coefficient and the 
surface wind. 

5.2 Heat Sources/Sinks 

Sensible heating is evaluated for the lowest model 
level over both ocean and land areas. Over the sea surface, 
the sensible heat flux is evaluated by the bulk method and 
the coefficient of sensible heat flux is treated as a func- 
tion of the static stability and wind speed. A Bowen ratio 
(which depends on the temperature and relative humidity of 
the lowest layer) relationship is used over land areas. 
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A moist convective adjustment scheme is applied at grid 
points where the relative humidity is greater than 70% 
in order to estimate the heating due to the release of 
latent heat and the convective transport of heat. A dry 
convective adjustment procedure is applied at the remaining 
grid points to maintain hydrostatic stability. 

5.3 Moisture Sources/Sinks 

Evaporation from the surface (ocean and land) into the 
lowest model layer is treated m a manner similar to the 
sensible heat calculation described in Section 5.2. 

A moist convective adjustment scheme is applied at 
points where the relative humidity is greater than 70% m 
order to estimate the precipitation and the convective 
transport of water vapor. 

6 . 0 Inputs 

The forecast model initial data is obtained from 
analyses of the mass structure, the moisture, the surface 
pressure and the sea surface temperature. 

Initial winds are obtained by combining the non- 
divergent component obtained from the stream function through 
the balance equation using the geopotential , and the 
irrotational component obtained from the velocity potential 
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as calculated from the quasi-geostrophic vertical p-velocxty 
satisfying the continuity equation. 

7 . 0 Outputs 

The twenty-four hour prognostic charts produced from 
the FLM forecast include: 

a) surface pressure, 

b) 24-hour rainfall, 

c) dew-point depression at 700 mb, 

d) vertical p-velocity (omega) at 700 mb, 

e) relative vorticity (zeta) at 700 mb, and 

f) geopotential contours and isotherms at 850 mb 
and at 500 mb. 

8 . 0 Computer Information. 

The 6L-FLM forecast model is run on a HITAC 8700/8800 
Computer System and requires about 27-28 minutes to make a 
one-day forecast. 

9.0 References 
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